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Abstract 


This  report*  summarizes  the  3  year  experimental  and  theoretical  research  program  conducted  by 
the  Fracture  Research  Laboratory  under  AFOSR  support  and  devoted  to  some  aspects  of  brittle 
fracture**. 

The  behavior  of  materials  subjected  to  brittle  fracture  is  studied  using  numerical  simulations  of 
subcritical  crack  growth.  Various  types  of  loadings  are  considered,  i.e.,  simple  and  eccentric 
tension,  three  and  four  point  bending,  dipole  on  the  crack  boundaries,  constant  and  time  varying 
loads.  As  a  result,  lifetime-load  relations  are  established,  a  new  approach  to  brittle  material 
toughness  assessment  is  proposed,  and  a  method  of  reliability  prediction  for  stmcmral 
components  is  developed. 


*  This  report  is  the  Ph.D.  thesis  of  Daniel  Baron,  which  was  defended  in  August  1996. 

**  The  other  works  in  the  framework  of  this  program  were  previously  reported: 

L  Toughening  of  Brittle  Materials,  1994. 

2.  The  Time  Dependency  of  the  Necking  Process  in  PC,  1995. 
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SUMMARY 


In  engineering  applications,  plastic  structural  components  are  often  needed  to  perform  for 
many  years  under  relatively  small  stresses.  These  components  usually  fail  suddenly  and  without 
warning,  when  a  slowly  growing  crack  originating  from  a  material  defect  reaches  a  critical  length, 
and  then  accelerates,  causing  almost  instantaneous  fracture.  In  order  to  determine  if  a  certain 
plastic  is  adequate  for  a  particular  structural  application,  it  would  be  ideally  desirable  to  make 
many  lifetime  tests  of  the  plastic  component  functioning  in  the  actual  mechanism,  and  then  to 
build  a  probability  distribution  of  component  lifetime.  The  distribution  would  provide  a  measure 
of  the  plastic’s  adequacy.  This  method  is  in  effect  used,  when  the  failures  of  a  commercially 
released  mechanism  are  tracked.  But  the  knowledge  comes  with  the  cost  of  real  world  failures.  It 
is  seldom  if  ever  possible  to  conduct  failure  tests  of  a  mechanism  which  is  being  designed.  For 
instance,  if  a  new  machine  was  being  designed  which  required  a  plastic  component,  failure  testing 
of  the  component  would  require  that  many  of  the  machines  be  built  and  each  operated  until  its 
component  failed.  The  time  to  conduct  the  failure  tests  would  be  comparable  to  the  machine’s 
design-lifetime,  and  so  for  a  long  term  design-lifetime,  completing  the  tests  would  be  greatly 
impractical. 

The  usual  alternative  to  mechanism  design-lifetime  failure  tests,  is  to  perform  failure  tests 
on  plastic  samples  using  material  testing  machines.  The  use  of  testing  machines  eliminates  the 
need  to  use  the  actual  mechanism  during  the  tests.  Standard  tests  are  performed  in  order  to 
provide  a  way  to  compare  plastics,  and  also  to  estimate  lifetimes.  The  standard  tests  are  performed 
using  standard  specimen  geometries.  The  actual  geometry  of  the  component  which  is  being 
designed  is  not  accounted  for.  One  or  more  design  parameters  (initial  crack  length,  applied  stress, 
temperature)  is  increased  enough  above  its  design  value,  so  that  the  testing  can  be  completed 
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within  the  allotted  time.  Empirical  methods  are  used  to  extrapolate  from  the  determined  test 
lifetimes  to  (much)  longer  estimated  design  hfetimes. 

The  Crack  Layer  Model  (CLM)  for  slow  crack  growth  in  engineering  plastics  was  created 
by  Professor  A.  Chudnovsky  in  1984  [1],  The  model  provides  a  mechanical  explanation  of  slow 
crack  growth.  A  numerical  implementation  of  the  equations  of  the  model  can  be  used  to  predict 
the  time  evolutions  of  a  crack  and  the  damage  zone  which  precedes  it,  in  a  plastic  structural 
component.  Computer  modeling  of  a  particular  problem  requires  the  specification  of  certain  of  the 
plastic’s  material  properties.  The  values  of  these  properties  are  found  by  performing  standard 
short  term  tests  on  material  specimens.  Also  required  is  the  ability  to  calculate  the  stress  intensity 
factor  (SEE)  at  the  crack  tip  for  the  specific  component  geometry  and  loading.  If  the  required 
material  properties  and  SIF  calculation  are  available,  then  the  numerical  solution  is  found  by 
solving  a  system  of  two  ordinary  differential  equations  thru  time.  A  lifetime  for  the  component  is 
found  without  resorting  to  historical  data,  long  term  mechanism  failure  tests,  or  short  term 
material  tests  which  use  empirical  manipulations. 

This  thesis  is  basically  a  numerical  implementation  of  the  equations  of  the  CLM. 
Algorithms  are  presented  for  the  cases  of  time  varying  prescribed  load  and  time  varying 
prescribed  displacement.  Simplified  algorithms  are  also  presented,  which  save  much  computing 
energy  and  under  certain  conditions  give  approximately  the  same  answers.  All  of  the  algorithms 
apply  to  the  simple  circumstance  of  a  plane  single  edge  notched  (SEN)  specimen,  with  a  straight 
crack.  A  summary  of  the  theory  of  the  CLM  is  also  presented.  The  output  from  computer 
simulations  employing  the  presented  algorithms  are  reported  and  analyzed.  Computed  types  of 
crack  growth  are  compared  to  experimentally  observed  types.  Computed  applied  stress  -  lifetime 


SUMMARY  (continued) 


relations  are  compared  to  experimentally  found  relations.  Computed  applied  stress  -  lifetime 
relations  are  compared  for  different  types  of  loadings.  A  method  is  presented  for  calculating  the 
reliability  function  (a  function  of  time  which  gives  the  probability  that  failure  has  not  happened) 
of  a  plastic  structural  component.  A  computer  simulation  is  made  of  the  ASTM  method  for 
measuring  the  fracture  toughness  of  a  plastic.  It  is  shown  that  the  measurement  found  is 
displacement  rate  dependent,  and  because  of  this  the  proposal  is  made  that  fracture  toughness  is 
not  a  material  property. 
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INTRODUCTION 


This  thesis  is  another  in  a  sequence  of  investigations  by  students  of  Professor  A.  Chudnovsky 
into  the  behavior  of  the  Crack  Layer  Model  (CLM)  for  engineering  plastics.  Previous  theses 
include  those  of  W.L.  Huang  and  K.  Kadota.  The  research  for  this  thesis  consisted  of  the 
numerical  implementation  of  the  CLM  for  various  loading  types  (tension,  three  point  bending, 
etc.)  and  histories  (fixed  load,  time-varying  load,  time-varying  displacement).  The  algorithms 
which  were  used  to  construct  a  computer  program  are  specified.  The  output  from  the  computer 
simulations  which  were  made  are  described  and  analyzed.  A  method  is  shown  for  determining  the 
reliability  of  a  plastic  structural  part  for  any  value  of  service  time.  It  is  demonstrated  that  the 
accepted  method  for  calculating  material  fracture  toughness  is  rate  dependent,  and  therefore 
ambiguous.  A  proposal  is  made  that  fracture  toughness  is  not  a  material  property,  and  that  one 
material  can  be  found  to  be  more  fracture  resistant  than  another,  only  by  accounting  for  all  of  the 
parameters  which  comprise  a  particular  usage  environment. 

An  understanding  of  the  slow  crack  growth  (subcritical  crack  growth)  in  plastics  which  leads 
to  long-term  brittle  fracture  is  very  important  to  industry.  The  use  of  plastics  as  structural 
components  is  always  increasing.  A  crack  which  causes  failure  usually  grows  from  a  defect  which 
is  originally  present  in  the  material.  The  defect  (often  a  void)  acts  as  an  initial  crack  with  length 
equal  to  the  defect’s  largest  dimension.  When  subjected  to  stress,  the  crack  slowly  grows  from  this 
beginning  length.  For  temperatures  in  the  range  of  room  temperature,  slow  crack  growth  occurs  at 
relatively  low  stresses,  and  is  the  main  cause  of  failure  for  plastic  parts  which  are  in  service  for 
many  years.  The  time  period  from  the  beginning  of  service  until  failure  is  called  the  lifetime  (or 
time  to  failure)  and  denoted  by  r^.  An  understanding  of  slow  crack  growth  would  allow  the 
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development  of  an  accelerated  material  testing  procedure,  based  on  the  experimental 
determination  of  material  parameters  and  the  numerical  modeling  of  the  crack  growth. 

Experimental  data  of  long-term  fracture  due  to  both  constant  and  cyclic  loading  have  been 
reported  for  many  plastics  [2-6],  and  are  very  detailed  for  polyethylenes  [7-24].  For 
polyethylenes,  a  region  of  damaged  material  referred  to  as  the  process  zone  (PZ),  is  always 
present  in  front  of  the  crack  tip.  The  PZ  forms  almost  immediately  after  loading,  because  the 
crack  tip  singularity  magnifies  the  transverse  tensile  stress  in  the  region  so  that  the  material’s  yield 
(drawing)  stress  is  reached.  In  polyethylene  and  some  other  plastics  (for  instance,  polymethyl 
methacrylate)  microvoids  develop  and  cause  the  PZ  to  consist  of  disconnected  drawn  fibers  which 
are  perpendicular  to  the  direction  of  the  crack.  The  behavior  of  PZ  material  is  usually  considered 
to  be  similar  to  that  of  drawn  (necked)  material  [25,26].  It  has  been  shown  that  the  time  dependent 
properties  of  the  PZ  fibers  (i.e.,  disentanglement  and  creep)  greatly  influence  the  crack  growth 
rate  and  lifetime.  Analyzing  the  PZ  from  the  viewpoint  of  continuum  mechanics  is  usually  based 
on  the  Dugdale-Barenblatt  (DB)  model  [27,28].  However,  observations  have  shown  that  actual 
PZs  have  significantly  shorter  lengths  than  are  predicted  by  the  DB  model  [29]. 

Experiments  have  demonstrated  that  with  respect  to  time,  slow  crack  growth  may  be 
continuous  (monotonically  increasing)  or  discontinuous  (intermittent).  Observations  of 
discontinuous  growth  in  polyethylenes  subjected  to  fatigue  loading  have  existed  for  a  long  time 
[3],  but  only  recently  has  discontinuous  growth  been  observed  for  constant  loading  [19]. 
According  to  [19],  the  growth  is  discontinuous  because  the  crack  is  arrested  until  the  long  chain 
molecules  of  the  drawn  PZ  fiber  in  front  of  the  crack  tip  untangle  and  then  the  fiber  breaks.  The 
rate  of  disentanglement  for  a  drawn  fiber  is  a  material  property. 

Descriptions  of  the  experimentally  observed  crack  kinetics  have  been  based  on  the  assumption 
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that  the  stresses  and  strains  in  the  neighborhood  of  the  crack  tip  completely  determine  the 

conditions  for  crack  initiation  and  growth.  This  means  that  the  rate  of  crack  growth  (/)  (which  is 
sometimes  approximately  measured  experimentally  by  using  the  rate  of  the  notch  tip  opening 

displacement  (t{N) )  [18,19])  is  a  function  of  only  such  fracture  parameters  as  the  stress  intensity 
factor  (K),  /-integral,  etc.  Observations  [7-22]  have  allowed  the  approximate  determination  of 
the  crack  growth  rate  under  constant  load  as  linearly  proportional  to  a  power  of  the  stress  intensity 

factor  at  the  crack  tip,  i.e.,  /  =  ,  t{N)  =  C2K\  where  Cj ,  C2  and  j  are  material  constants. 

According  to  [30],  for  polyethylenes  j  varies  from  3  to  5.  The  empirical  relations  mentioned  have 
been  established  for  the  simplest  conditions  of  crack  propagation,  i.e.,  the  specimen  material  is 
homogeneous  and  the  crack  path  is  straight.  Study  of  the  case  when  the  specimen  contains  an 
inclusion  (a  region  with  different  material  properties)  has  shown  that  crack  kinetics  cannot  be 
completely  determined  from  the  crack  tip  parameters.  Additional  parameters  which  account  for 
the  effect  of  the  PZ  on  the  crack  direction  and  growth  rate  must  be  specified. 

With  respect  to  lifetime  prediction,  the  most  important  experimental  result  has  been  the 
development  of  relations  between  lifetime,  applied  stress  and  temperature.  Data  about  the 
lifetimes  for  internally  pressurized  high  density  polyethylene  pipes  at  various  load  levels  and 
temperatures  are  reported  in  [5].  Experiments  [16-18]  on  single-edge  notched  (SEN)  tensile 
specimens  of  polyethylenes  with  various  molecular  weight  distributions  and  branch  densities  have 
allowed  the  determination  of  lifetime  as  a  function  of  stress  level,  temperature  and  notch  length. 
For  a  constant  load  and  a  particular  material,  specimen  geometry,  load  type  (tension,  three-point 

bending,  etc.)  and  temperature,  lifetime  is  a  function  of  the  stress  level  a  only,  i.e.,  a  ” .  In 
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[5]  a  is  the  pipe  hoop  stress,  and  in  [16-18]  it  is  the  applied  tensile  stress,  a  is  a  material  constant 
which  depends  little  or  none  on  temperature,  and  varies  from  approximately  2.5  to  5.5  for 
different  polyethylenes.  It  has  been  shown  that  a  particular  treatment  of  the  lifetime-stress  relation 
obtained  from  short-term  fracture  tests  (at  high  temperatures)  can  give  a  reasonable  estimate  of 
the  lifetime-stress  relation  for  long-term  processes  at  room  temperature  [16-18,23,24]. 

There  are  now  three  approaches  used  for  modeling  slow  crack  growth  in  plastics.  The  first 
assumes  that  the  kinetics  are  completely  determined  by  the  time-dependent  properties  of  the  bulk 
material  which  surrounds  the  crack  and  PZ.  According  to  the  most  detailed  theory  of  this  kind 
[31],  the  crack  grows  because  of  viscoelastic  deformation  of  the  bulk  material.  The  second 
approach  connects  crack  growth  to  the  time-dependent  processes  which  take  place  in  the  PZ 
material.  For  polyethylenes,  the  time-dependent  processes  are  disentanglement  and  creep  of  the 
PZ  fibers.  This  approach  has  been  studied  experimentally  in  [32-42].  An  explanation  of  notch  tip 
opening  displacement  as  being  due  to  visco-elongation  of  the  PZ  fibers  is  given  in  [30].  The  third 
approach  (which  this  thesis  follows)  considers  the  PZ  region  to  consist  of  damaged  material,  and 
the  crack  and  PZ  as  the  two  co-dependent  components  of  a  system  named  the  crack  layer  (CL) 
[1,43].  CL  kinetics  are  analyzed  using  irreversible  thermodynamics,  and  a  complete  solution  for 
the  general  problem  has  been  found.  An  approximation  of  the  CLM  has  permitted  numerical 
modeling  of  various  slow  crack  growth  phenomena,  including  discontinuous  crack  growth  in 
polyethylenes  [29,44-45],  etc.  The  approximation  used,  is  to  assume  the  CL  to  be  a  two  parameter 
system,  the  crack  length  (/)  and  the  crack  layer  length  (L)  (L  is  equal  to  the  sum  of  /  and  the 
active  PZ  length  (Z^)).  (The  numerical  algorithms  which  are  used  for  this  thesis,  use  the  two 


parameter  approximation.) 
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The  approximation  of  the  CLM  which  was  just  stated,  is  plausible  because  observations  of 
actual  process  zones  in  plastics  have  shown  that  the  process  zone  length  (Z^)  is  always  much 

greater  than  the  process  zone  width  Therefore,  in  the  equations  of  the  thermodynamic 
analysis,  the  terms  containing  are  relatively  very  small,  and  when  the  approximation  is  used, 

are  neglected.  The  bulk  material  is  assumed  to  be  linearly  elastic.  (Actually,  all  engineering 
plastics  possess  viscoelastic  properties,  but  the  criteria  for  crack  and  PZ  growth  are  related  only  to 
the  elastic  energy  of  the  material.  For  the  case  of  prescribed  load,  viscoelastic  deformation  does 
not  affect  the  stress  field  and  so  does  not  affect  lifetime.  For  the  case  of  prescribed  displacement, 
viscoelastic  deformation  would  affect  the  stress  field,  and  so  lifetime  would  be  affected.  However, 
viscoelasticity  is  outside  the  boundary  of  this  thesis.)  During  CL  growth,  the  potential  energy  of 
the  bulk  material  is  used  to  break  the  process  zone  fibers  and  extend  the  crack.  As  the  crack 
extends,  its  stress  intensity  factor  increases,  which  forces  the  PZ  to  extend.  The  time-dependent 
changes  of  the  PZ  region  are  collectively  called  material  degradation.  For  prescribed  external 
forces  and  fixed  temperature,  the  Gibbs  free  energy  (G)  of  the  bulk  material  and  PZ  is  the 
potential  which  indicates  deviation  of  the  system  from  equilibrium.  From  irreversible 
thermodynamics,  the  thermodynamic  force  for  crack  length  extension  (X^)  is  defined  as  the 

partial  derivative  of  G  with  respect  to  the  crack  length  (/),  and  the  thermodynamic  force  for  crack 
layer  length  extension  )  is  defined  as  the  partial  derivative  of  G  with  respect  to  the  crack  layer 

length  (L).  A  crack  layer  parameter  (/,  L)  grows  whenever  the  value  of  its  thermodynamic  force 
is  positive.  Numerical  solution  of  the  CLM  kinetic  equations  has  demonstrated  emulation  of  both 
types  of  experimentally  observed  slow  crack  growth,  continuous  and  discontinuous.  Stability 


analysis  of  the  CLM  indicates  that  growth  of  /  and/or  L  occurs  whenever  the  system  is  in  an 
unstable  state.  Growth  continues  until  a  new  stable  state  is  attained.  For  a  discontinuous  process, 
stable  states  exist,  and  growth  stops  whenever  one  is  reached.  The  time-dependent  degradation  of 
the  PZ  material  causes  a  stable  system  to  become  unstable.  For  a  continuous  process,  there  are  no 
possible  stable  states,  so  growth  is  never  suspended. 

This  thesis  is  about  the  numerical  analysis  of  the  CLM  for  various  types  of  loads  and  load 
histories.  The  CLM  was  numerically  implemented  by  developing  a  Fortran  computer  program. 
Chapter  1  contains  a  brief  review  of  the  CLM  theory,  and  a  derivation  of  the  kinetic  equations 
which  are  used  to  model  slow  crack  growth  in  engineering  plastics.  Chapter  2  relates  the  detailed 
numerical  algorithms  which  were  used  to  construct  the  program  which  solves  the  equations  thru 
time.  Algorithms  are  given  for  the  cases  of  prescribed  fixed  load,  prescribed  time-varying  load  and 
prescribed  time-varying  displacement.  Simplified  algorithms  are  also  given,  which  are  usually 
very  accurate  when  modeling  discontinuous  crack  growth.  The  simplified  algorithms  are  included 
because  their  use  can  greatly  reduce  the  work  which  the  computer  must  do.  Chapter  3  analyzes  the 
output  from  the  fixed  load  computer  simulations  which  were  made.  The  thermodynamic  forces 
are  separated  into  forcing  and  resisting  parts.  The  set  of  simulation  input  parameters  is  presented 
in  dimensionless  form,  and  the  requirements  for  simulation  similarity  are  indicated.  The  factors 
which  determine  whether  a  simulation  is  continuous  or  discontinuous  are  discussed.  The  output 
from  the  simplified  algorithm  is  compared  to  the  output  from  the  complete  algorithm.  Lifetime- 
applied  load  relations  are  provided  for  both  continuous  and  discontinuous  simulations,  and  also 
for  various  types  of  loading.  Chapter  4  shows  a  way  to  calculate  the  reliability  function  for  a 
plastic  structural  component.  The  method  involves  performing  destructive  tests  on  material 
specimens  in  order  to  build  a  distribution  of  the  largest  defects,  and  then  performing  computer 
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simulations  to  develop  a  lifetime-defect  size  relationship.  Chapter  5  uses  simulations  of 
prescribed  time-varying  displacement  to  imitate  the  ASTM  method  for  measuring  the  fracture 
toughness  of  a  material.  It  is  shown  that  the  results  depend  on  displacement  rate,  and  so  can  be 


misleading. 


1.  THE  CRACK  LAYER  MODEL  FOR  SLOW  CRACK 
GROWTH  IN  ENGINEERING  PLASTICS 

1.1  Introduction 

In  this  thesis,  the  modeling  of  slow  crack  growth  which  causes  brittle  fracture  of 
engineering  plastics,  is  based  on  the  idea  that  slow  crack  growth  is  the  evolution  of  two  elements. 
The  first  element  is  the  crack,  and  the  second  is  a  region  of  damaged  material  preceding  the  crack 
tip  which  is  called  the  process  zone  (PZ).  Following  [1],  the  system  consisting  of  the  crack  and  PZ 
is  called  the  crack  layer  (CL).  Advance  of  the  CL  results  from  the  interactions  between  the  crack 
and  the  PZ,  and  between  the  PZ  and  the  body  surrounding  it,  which  is  called  the  matrix.  The 
equations  governing  CL  growth  are  derived  using  the  conventional  procedure  of  the 
thermodynamics  of  irreversible  processes.  The  temperature  and  load  are  considered  to  be 
prescribed,  so  the  relevant  thermodynamic  potential  is  Gibbs  free  energy.  Gibbs  free  energy 
accounts  for  the  potential  energy  of  the  matrix  and  mechanism  maintaining  the  load,  the  energy 
required  for  PZ  formation,  and  the  energy  accumulation  on  new  crack  boundaries.  Partial 
derivatives  of  Gibbs  potential  with  respect  to  the  craek  length  and  PZ  length  give  the 
thermodynamic  forces  for  crack  and  PZ  extensions.  Kinetic  equations  are  formulated  as  linear 
relations  between  the  thermodynamic  forces  and  the  rates  of  crack  length  extension  and  PZ  length 
extension.  The  process  becomes  unstable  when  the  rates  begin  to  monotonically  increase.  The 
time  period  from  load  application  until  CL  instability  is  called  the  lifetime  or  time  to  failure,  and 
written 

Only  the  opening  mode  of  crack  loading  (mode  1)  and  a  straight  crack  are  considered.  In  this 
case  the  PZ  of  an  engineering  plastic  has  the  shape  of  a  thin  strip  ahead  of  the  crack  tip  and  along 
the  crack  centerline.  Therefore,  the  CL  can  be  described  by  just  two  independent  parameters,  the 
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crack  and  PZ  lengths.  A  system  of  two  kinetic  equations  is  formed.  The  system  is  nonlinear 
beeause  the  thermodynamic  forces  depend  nonlinearly  on  the  two  parameters,  and  because 
constraints  on  the  rates  of  crack  and  PZ  extension  are  expressed  as  inequalities. 

1.2  Energy  balance  for  crack  layer  growth.  Thermodynamic  forces 

It  is  usually  accepted  that  the  PZ  material  is  similar  to  the  material  obtained  from  cold 
drawing  (necking)  the  original  plastic.  So,  the  PZ  material  is  supposed  to  be  formed  under  the 
drawing  stress  and  stretched  according  to  the  draw  ratio  (X).  and  X  are  material 

properties  of  a  particular  plastic. 

What  X  means  will  now  be  explained.  If  a  plastic  specimen  is  loaded  in  tension  along  the 
longitudinal  axis,  when  is  reached  a  transverse  band  of  necked  material  will  appear  at  the 

location  of  the  specimen’s  weakest  cross  section.  The  band  thickness  will  grow  from  zero.  While 
the  load  is  maintained,  the  band  thickness  will  increase  as  its  boundaries  (the  two  end  cross 
sections  of  the  band)  advance  in  both  longitudinal  directions.  The  boundaries  can  be  thought  of  as 
material  transformers.  As  a  boundary  moves  across  a  strip  of  unnecked  material,  the  strip  is 
transformed  into  necked  material.  The  axial  strain  in  the  unnecked  material  will  be  called  ,  and 

the  axial  strain  in  the  necked  material  will  be  called  Eq.  It  is  important  to  note  that  Eq  is  a 
constant,  and  that  once  necking  begins  E^  becomes  constant.  (It  should  also  be  noted  that  Eq  is  a 
constant  only  until  all  of  the  material  of  the  specimen  becomes  necked.  When  no  unnecked 
material  remains,  and  if  the  load  is  still  maintained,  Eq  will  begin  to  grow.)  An  arbitrary 

longitudinal  interval  measured  on  the  unloaded  specimen,  will  have  length  /j .  When  the  specimen 
is  loaded  until  necking  occurs,  the  interval  will  have  a  new  length.  If  the  measured  interval  is 
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completely  outside  the  band  of  material  which  becomes  necked,  /j  will  become  I2.  If  the 
measured  interval  is  completely  inside  the  band  of  material  which  becomes  necked,  /j  will 
become  13 .  X,  is  defined  as  the  ratio  of  to  I2  ■ 


Therefore, 


So, 

Eq  =  X(l  +  Ey)  -  1  . 

Since, 

Eu  «  1  , 

Gq  ~  ^  ^  • 

The  approximation  of  equation  1.5  will  be  used  in  the  remainder  of  this  section,  i.e., 
Eq  =  X  -  1 . 


(1.2) 

(1.3) 

(1.4) 

(1.5) 

(1.6) 


Since,  as  was  stated  in  section  1.1,  the  PZ  is  shaped  like  a  thin  strip,  the  stress-strain  state  of 
the  matrix  due  to  the  presence  of  the  CL  can  be  approximated  by  a  slit  (a  straight  one  dimensional 


11 


cut  in  the  matrix)  with  tractions  of  magnitude  distributed  along  the  location  of  the  PZ  (Figure 

1). 
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Figure  1.  One  dimensional  approximation  of  crack  layer 


The  energy  balance  for  PZ  formation  will  be  derived.  Gibbs  potential  includes  the  change  in 
the  potential  energy  of  the  matrix  and  mechanism  maintaining  the  prescribed  load,  and  the  change 
in  the  deformation  potential  of  the  material  which  becomes  the  PZ.  The  change  in  the  total 
potential  energy  of  the  matrix  and  load  is 
U=F-  W, 

m  m  < 


(1.7) 
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is  the  change  in  matrix  elastic  energy,  and  is  the  work  done  by  the  applied  load  moving 
thru  the  displacement  caused  by  PZ  formation.  According  to  numerous  experimental 
observations,  cold  drawing  (necking)  of  plastics  gives  the  stress-strain  function  (a(8))  shown  in 
Figure  2.  The  function  is  considered  to  be  a  material  property  for  the  PZ  during  its  formation.  The 
change  in  the  PZ  deformation  potential  is  defined  as 

Fpz=  J  hix)\  J  o(e)JeWx.  (1.8) 

x=l  ^  e=0  ^ 


h{x)  is  the  width  of  the  original  material  strip  which  is  transformed  into  the  PZ.  The  crack  layer 
opening  displacement  (CLOD)  is  defined  as 

6p^(^)  =  e(x)/i(x).  (1.9) 

From  equation  1.6, 


h{x)  = 


^pz(x) 

x-1  • 


(1.10) 


Define 


e  =  X-l 

Q=  J  CT(£)j£.  (1.11) 

e  =  0 

Q  is  the  area  under  a  =  a(£)  in  Figure  2.  Equation  1.8  can  be  written  as 

x  =  L 

=  ^  I  (1.12) 

x  =  l 

Define 


a'^G^,(X-i). 


(1.13) 
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Q'  is  the  area  under  a  =  in  Figure  2.  If  the  PZ  material  was  perfectly  plastic,  then  equation 
1.12  would  become 


X 

q: 

X-  1 

JC 


=  L 

I  dp^{x)dx. 
=  l 


(1.14) 


Figure  2.  Process  zone  formation  stress-strain  diagram 


(The  definition  of  £2’  implies  that  Young’s  modulus  is  infinite,  but  the  resulting  error  in  is 
very  small.)  Define 
-  Q-Q’ 


(>0) 


(1.15) 
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y  is  the  ratio  of  the  overloading  energy  to  the  perfectly  plastic  energy  during  PZ  formation.  From 
equations  1.12-1.15, 

Fpz  =  (1-16) 

Define 


n„'  =  n  +  F  . 

m  m  pz 

The  change  in  Gibbs  potential  is 

From  applying  Clapeyron’s  theorem  (Figure  1), 

{X  =  L  X  =  L 

j  J  <5^{x)h{x)dx 

x=l  x=0 

and 


(1.17) 


(1.18) 


(1.19) 


(1.20) 


(Within  each  integral  of  equation  1.19,  5(x)  is  the  displacement  of  the  particular  traction.) 
The  differential  of  Gibbs  potential  due  to  CL  growth  is 


-dG  = 


bl  bl 


dL  ,  (>0). 


(1.21) 


Equation  1.21  is  equal  to  zero  when  the  CL  is  in  equilibrium.  2y  is  the  energy  required  for  the 
crack  to  extend  into  the  PZ,  i.e.,  the  Griffith  fracture  energy  of  the  damaged  material.  (Note  that  y 
and  2y  are  not  related.)  Define 


an^'  BF" 
^  ^ _ ^ _ 


Bl 


Bl 


-2y, 


(1.22) 


and 


Energy  balance  1.21  is  written 


-dG  =  Xidl  +  X^dL  ,  (>0).  (1.24) 

X[  is  called  the  driving  force  (thermodynamic  force)  for  crack  length  extension,  and  X^  is  called 

the  driving  force  (thermodynamic  force)  for  crack  layer  length  extension.  By  using  the  solution  of 
the  elastic  problem  shown  in  Figure  1,  the  partial  derivatives  in  equations  1.22  and  1.23  can  be 
expressed  in  terms  of  conventional  fracture  mechanics  parameters. 


■  3/  -  +  5^r) 

(1.25) 

3n„,’  1  2 

■  3L  = 

(1.26) 

(1.27) 

dF"  2yK^^ 

BL  =  E 

(1.28) 

E  is  Young’s  modulus.  6^  is  the  crack  opening  displacement  (COD)  at  the  crack  tip  (x  =  /),  due 
to  the  influence  of  a^(j:)  .  5^^  is  the  COD  at  =  /,  due  to  the  influence  of  is  the 

stress  intensity  factor  (SIF)  at  x  =  L  due  to  the  influence  of  a^{x) .  is  the  SIF  at  x  =  L 
due  to  the  influence  of  .  Equations  1.22  and  1.23  are  written 

Xi  =  -cy^^(8^  +  (l +y)8^,)-2Y,  (1.29) 

+  ( 1  +  2y)^^^)  .  (1.30) 
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So,  the  driving  forces  are  expressed  in  terms  of  the  standard  parameters  of  linear  fracture 
mechanics. 

Each  of  the  driving  forces  can  be  separated  into  two  parts.  The  forcing  parts  of  the  driving 
forces  are  defined  as 


(1.31) 

(1.32) 

The  resisting  parts  of  the  driving  forces 

are  defined  as 

(1.33) 

(1.34) 

Each  of  the  quantities  Dj,  D^,  Ri, 

are  always  positive.  Equations  1.29  and  1.30  are  written 

Xi  =  Di-Ri, 

(1.35) 

(1.36) 

The  crack  length  will  increase  only  when  X^  is  positive.  The  crack  layer  length  will  increase  only 

when  is  positive. 

1.3  DF.nRAPATION  OF  THE  PROCESS  ZONE  MATERIAL 

The  two  parameter  model  [38-40]  assumes  that  resistance  of  the  PZ  material  to  crack 
extension  decreases  with  time.  The  PZ  material  is  said  to  degrade  with  time.  PZ  formation  occurs 
under  the  drawing  stress  and  with  stretching  subject  to  the  draw  ratio  (X,).  The  PZ  material 

is  incompressible,  so  the  true  stress  acting  on  the  PZ  is  defined  to  be  .  The  surface  energy 
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(2y)  represents  a  measure  of  the  PZ  resistance  to  crack  extension.  The  PZ  consists  of  many  fibers 
perpendicular  to  the  direction  of  the  crack  which  are  supposed  to  be  drawn. 

If  a  PZ  fiber  is  drawn  at  time  t  =  0,  it  has  an  initial  axial  drawing  strain  (Eq).  As  time  (0 

increases,  and  the  fiber  continues  to  be  subjected  to  stress  the  fiber  axial  strain  (e)  will 

increase  until  the  critical  fiber  axial  strain  (e^)  is  reached,  and  the  fiber  breaks.  (The  crack 
advances  whenever  the  fiber  at  the  crack  tip  breaks.)  The  critical  axial  strain  (£^ )  is  reached  at  the 
critical  time  of  drawing  and  are  material  constants.  The  axial  strain  (e)  grows  from  Eq 

to  E^  according  to  a  material  function  e(t)  ■  Define 

AE(0  =  E(0-eo,  (1-37) 
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A  PZ  fiber  strain-time  relation  is  shown  in  Figure  3. 


Figure  3.  Process  zone  fiber  strain-time  relation 


At  time  t  =  0  the  work  done  on  the  PZ  fiber  is 

At  time  t  the  work  done  on  the  PZ  fiber  is 
w(t)  =  WQ  +  X<T^^AE(t). 

At  rupture  time  t  =  the  amount  of  work  done  on  the  PZ  fiber  is 


(1.39) 


(1.40) 


(1.41) 
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The  surface  energy  required  for  crack  length  extension  ,  see  equation  1.29)  is  defined  as 

ly{t)  =  w^-w{t) .  (1-42) 

2y{t)  is  the  amount  of  energy  required  to  break  the  PZ  fiber  at  time  t .  The  initial  and  maximum 
value  of  the  surface  energy  (2yQ )  is  a  material  constant  and  defined  as 

2Yo  =  2y(0)  =  (1-43) 

2y{t)  can  be  written 

2y«)  =  2Yo(i-^).  (1-44) 

If,  for  instance  the  material  function  Ae(r)  is  linear,  i.e.,  if 

Ae(0  =  at,  (1.45) 

then 

If  Ae(0  is  exponential,  i.e.,  if 

A£(0  =  Ae,(l-e"P'),  (1.47) 

then 

2y(0  =  2Yoe"P'  ,  (0<r<oo).  (1.48) 

In  equation  1.48,  t  represents  the  elapsed  time  since  a  PZ  fiber  was  drawn.  A  PZ  fiber  is  formed 
(and  drawn),  when  its  material  becomes  part  of  the  PZ.  So  later  in  this  thesis,  2y{Atp^)  will  be 

used  instead  of  2y{t).  signifies  the  elapsed  time  since  the  material  at  the  location  of  the 


crack  tip  became  part  of  the  PZ.  (The  interest  is  always  with  the  material  at  the  crack  tip,  because 
the  state  of  that  material  determines  whether  the  crack  can  grow,  i.e.,  in  order  for  the  crack  to 
grow,  the  PZ  fiber  directly  ahead  of  the  crack  tip  must  be  broken.) 

Note  that  another  material  constant  (^q)  will  be  used  later  in  this  thesis,  /cq  is  the  reciprocal 

of  the  time  until  rupture  (t^)  of  drawn  material,  i.e.. 


1.4  Kinetic  equations 


(1.49) 


The  relations  between  driving  forces  and  instantaneous  growth  rates  are  assumed  to  be 

linear. 


/  =  k^X„ 

(1.50) 

L  =  k^X^. 

(1.51) 

kj  and  ^2  are  material  constants  which  are  determined  experimentally.  /  is  the  instantaneous 
growth  rate  for  the  crack  length.  L  is  the  instantaneous  growth  rate  for  the  crack  layer  length  (for 

the  front  end  of  the  PZ).  If  a  calculated  value  of  /  or  L  is  ever  negative,  then  that  value  is  set  equal 
to  zero.  Doing  this  prevents  contradictions  with  physical  reality. 

Enforced  limitations  on  the  values  of  I  and  L  are 
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is  the  equilibrium  value  for  the  crack  layer  length.  For  a  particular  value  of  I ,  is  the  value 
of  L  which  causes  (equation  1 .30)  to  equal  zero.  Also,  for  a  particular  value  of  I ,  is  the 
maximum  possible  value  of  L . 

The  initial  conditions  for  the  system  of  equations  1.50-1.51  are 


/(r  =  0)  =  Zq, 

(1.53) 

o 

VI 

o 

II 

O 

II 

(1.54) 

The  two  degree  of  freedom  system  of  equations  1.50-1.51  can  be  solved  numerically  thru 
time,  yielding  the  evolutions  of  I  and  L . 

1.5  Conclusions 

The  driving  (thermodynamic)  forces  for  crack  extension  (X^)  and  crack  length  extension 
(X^)  were  derived.  The  opening  mode  (mode  1)  of  crack  loading  and  a  straight  crack  were 

assumed.  Therefore  it  was  guaranteed  that  the  length  of  the  PZ  would  be  much  larger  than  its 
width.  Having  a  narrow  PZ  allowed  the  CL  to  be  described  one  dimensionally,  and  by  just  two 
parameters,  the  crack  length  (  / )  and  the  crack  layer  length  (L ). 

The  definition  for  the  material  surface  energy  (or  degradation)  function  was  derived.  It  was 
assumed  that  the  PZ  was  composed  of  fibers  which  had  been  drawn.  This  allowed  the  strain  in  the 
fibers  of  a  newly  formed  PZ  to  be  calculated  using  the  material’s  draw  ratio  (A. ).  In  the  calculation 
of  the  area  under  a  =  in  Figure  2,  the  approximation  of  Young’s  modulus  {E)  being  infinite 
was  used. 

The  system  of  two  (nonlinear  ordinary  differential)  kinetic  equations  was  presented,  along 
with  the  initial  conditions,  and  the  constraints  on  the  calculated  values  of  the  crack  length  (/)  and 
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the  crack  layer  length  (L).  It  was  assumed  that  there  was  a  linear  relationship  between  a  driving 
force  and  the  associated  instantaneous  growth  rate. 


2.  NUMERICAL  ALGORITHMS  FOR  SOLUTIONS  OF 
CRACK  LAYER  MODEL  EQUATIONS 

2.1  Introduction 

This  chapter  will  show  numerical  procedures  for  determining  the  crack  length  l(t, ...), 
and  crack  layer  length  L(t, ...),  in  computer  simulations  of  specimen  lifetime  tests  using  the 
equations  of  the  Crack  Layer  Model  (CLM).  Algorithms  are  given  for  the  cases  of  fixed 
prescribed  load,  time-varying  prescribed  load,  and  time-varying  prescribed  displacement. 
Simplified  versions  of  the  algorithms  are  also  included.  They  are  based  on  the  implementation  of 
an  assumption  which  results  in  the  decoupling  of  equations  2.1.  The  associated  simplified 
simulations  require  little  computer  time  relative  to  simulations  using  the  complete  algorithms,  and 
give  approximately  the  same  output  for  the  modeling  of  discontinuous  crack  growth. 

The  goal  is  to  be  able  to  solve  the  system  of  two  nonlinear  ordinary  differential  equations 

i=hXi  (a)  ^2.1) 

L  =  k^Xi^  (b) 

i  is  the  instantaneous  rate  of  change  of  the  crack  length  with  respect  to  time.  L  is  the 
instantaneous  rate  of  change  of  the  crack  layer  length  with  respect  to  time.  A:j ,  a  material 

constant,  is  the  crack  length  kinetic  coefficient.  ^2 ,  a  material  constant,  is  the  crack  layer  length 
kinetic  coefficient.  X^  is  the  driving  force  for  crack  length  extension.  is  the  driving  force  for 
crack  layer  length  extension. 
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Xi=J,-2y 

{K^  +  Kj^)iK^  +  {l+2y)K,^) 


Xl  = 


(a) 

(b) 


(2.2) 


7j  is  the  rate  of  system  potential  energy  release  for  crack  length  extension.  2y ,  a  material 
function,  is  the  rate  of  surface  energy  required  for  crack  length  extension. 

Ji  =  ( 5^  +  ( 1  +  y )  5^,)  (2.3) 

6^  is  the  crack  tip  (at  a: -coordinate  x  =  1)  opening  displacement  due  to  the  remote  (applied) 
load,  y ,  a  material  constant,  is  the  specific  energy  of  drawing,  ,  a  material  constant,  is  the 
drawing  stress.  5^^  is  the  crack  tip  (at  x -coordinate  x  =  1)  opening  displacement  due  to  the 
process  zone  stress,  . 

2y  =  2yiAtp^)  (2.4) 

2y  is  a  decreasing  material  function,  of  the  elapsed  time  (At^^)  since  the  material  at  the  current 
location  (x-coordinate)  of  the  crack  tip  became  part  of  the  process  zone.  is  the  stress  intensity 
factor  at  L  (the  x-coordinate  of  the  front  of  the  crack  layer)  due  to  the  remote  load.  is  the 
stress  intensity  factor  at  L  due  to  the  process  zone  stress,  . 

2.2  Algorithm  for  SEN  specimen  sub.iected  to  creep  loading 

An  SEN  (single  edge  notched)  specimen,  is  a  plane  specimen  with  a  notch  cut  in  one  edge. 
The  notch  is  cut  perpendicular  to  the  specimen’s  longitudinal  axis.  “Creep  loading”,  means  the 
loads  to  which  the  specimen  is  subjected  are  constant. 
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2.2.1  Program  arrays 

The  computer  program  for  the  numerical  simulation  contains  arrays^^^  which  record 
the  values  of  the  main  variables  at  the  end  of  each  time  step.  These  arrays  include 


t[i] 

l[i]  ,  i  =  0,1,  (2.5) 

L[i] 

i  is  the  number  of  the  time  step,  and  n  is  the  total  number  of  time  steps  which  have  been 
completed,  in  order  to  reach  the  current  time,  r[n] .  The  initial  conditions  are  specified  as 


t[0]  =  0 

/[0]  =  /o  ,  (,Lo  =  Iq).  (2.6) 

L[0]  =  Lo 

2.2.2  Determination  of  the  remote  load  function 

The  remote  load  is  the  external  load  which  is  applied  to  the  specimen.  The  distance 
from  the  remote  load  to  the  crack  layer  is  assumed  to  be  large  enough,  such  that  St.  Venant’s 
principle  applies  with  respect  to  the  calculation  of  remote  load  stresses  in  the  region  of  the  crack 
layer.  The  remote  load  function,  o^(x) ,  is  the  normal  stress  which  would  be  present  in  the 


specimen  along  the  crack  centerline,  if  the  crack  layer  (the  crack  and  process  zone)  did  not  exist. 
For  instance,  for  a  plane  specimen  of  width  W  and  thickness  B ,  and  loaded  in  tension  by  force  F , 


a 


oo 


(2.7) 


For  a  plane  specimen  of  width  W,  thickness  B ,  and  loaded  in  the  plane  by  end  moments  M , 


(l)Note:  within  this  thesis,  a  symbol  immediately  followed  by  parentheses  “( )”  will  indicate  a  function,  and 
a  symbol  immediately  followed  by  square  braces  “[  ]”  will  indicate  an  array. 
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<Jco(^)  = 


(2.8) 


The  origin  of  the  x-axis  is  at  the  tension  edge  (which  is  also  the  notch  edge)  of  the  specimen. 

Note,  that  this  thesis  relates  only  to  the  states  of  plane  stress,  and  plane  strain,  so  that  B  is 
always  set  equal  to  one,  and  will  usually  not  be  mentioned. 

2.2.3  Initial  time  increment 

One  of  the  input  parameters  which  the  user  specifies  for  a  simulation,  is  the  initial  time 
increment  (A^q)  .  During  the  first  attempt  at  executing  every  time  step  (sometimes  an  attempt  will 


fail).  At ,  the  time  increment  used  during  the  time  step,  is  set  equal  to  Atp . 

Ar  =  A/q  (for  first  attempt  of  every  time  step)  (2.9) 

2.2.4  Calculation  loop  for  time  evolution  of  crack  and  crack  layer 

This  loop  is  repeated  for  each  time  step,  until  the  conditions  specified  within  the 
computer  program  which  indicate  failure  of  the  simulated  specimen  are  satisfied.  It  is  desired  to 
calculate  the  time  evolution  of  /(...),  the  x-coordinate  of  the  crack  tip,  and  of  L(...),  the  x- 
coordinate  of  the  front  of  the  crack  layer. 

2.2.4. 1  Increment  of  time  step  counter 

At  the  start  of  each  time  step,  the  time  step  counter  (n )  is  incremented, 

n  =  n+l.  (2.10) 

For  the  first  time  step. 


(n  =  n  +  1)  =>  (1  =  0  +  1). 


(2.11) 
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2.2.4.2  Determination  of  the  surface  enersx 

The  surface  energy  {2j(Atp^))  is  calculated  for  use  during  the  time  step.  The 

surface  energy  (Figure  4)  is  equal  to  the  amount  of  energy  required  to  extend  the  crack  a  distance 
such  that  one  area  unit  of  new  crack  surfaces  is  generated.  (For  a  specimen  with  unit  thickness 
{B  =  1 ),  this  distance  would  equal  one  half,  since  as  the  crack  extends  two  surfaces  are 
generated.)  When  the  crack  extends,  it  does  so  into  process  zone  material.  The  fibers  of  the 
process  zone  are  subjected  to  the  material  drawing  stress,  Subjection  to  causes  the 
material  of  the  fibers  to  degrade  and  the  energy  required  to  extend  the  crack  to  decrease,  as  time 
passes.  2y(  . . . )  is  a  monotonically  decreasing  material  function,  with  upper  limit  IJq  ,  and  lower 

limit  zero.  2Yo  is  equal  to  the  surface  energy  of  newly  drawn  and  still  undegraded  material. 

is  equal  to  the  elapsed  time  since  the  material  at  the  current  location  of  the  crack  tip  (at  / )  was 
drawn,  i.e.,  since  that  material  became  part  of  the  process  zone,  t  is  the  current  time,  t^^  is  the 

time  when  the  material  at  the  current  location  of  the  crack  tip  was  drawn,  i.e.,  when  that  material 
became  part  of  the  process  zone. 

=  (2-12) 
tp^  is  equal  to  the  time  which  makes  L{tp^, ...)  (the  x-coordinate  of  the  front  of  the  crack  layer  at 
time  tp^ )  equal  to  l(t, ...)  (the  x-coordinate  of  the  crack  tip  at  the  current  time  (t) ). 

L(tp^, ...)  =  lit, ...),  iv/iihtp^<t). 


(2.13) 
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In  order  to  determine  2y{Atp^) ,  it  is  necessary  to  find  .  The  current  time  t  is  known,  so  if 
can  be  found,  then  Ar„.  can  be  calculated.  In  order  to  determine  t  ,  it  is  first  necessary  to  find  j 

pz  pz 

such  that 

L[y]</[n-l]<L[;+l].  (2.14) 

(It  might  seem  strange  that  although  this  is  the  nth  time  step,  /(n  -  1)  is  used  in  equation  2.14. 
But  it  would  be  impossible  to  use  /(n) ,  since  it  has  not  been  calculated  yet.)  Finding  j  is  done  by 
searching  the  (ordered)  array  L[i] ,  beginning  with  L[0] ,  and  comparing  the  values  with 
/[n  -  1  ] .  Once  j  has  been  found,  then  by  linear  interpolation, 

^  ^ •  (2,15) 
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After  is  found,  can  be  calculated  using  equation  2.12.  Then  can  be  substituted  into 
the  prescribed  material  function  2 y(  . . . ) . 


a)  Exponential  degradation 

Figure  4.  Surface  energy  as  a  function  of  time 


2.2.4.3  Calculation  of  the  SIFs  and  CODs 

The  appendix  gives  the  algorithms  for  calculating  the  stress  intensity  factors 
(SIFs)  and  crack  opening  displacements  (CODs).  The  parameters  for  each  of  the  functions  will  be 


shown  here. 
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K„  =  KJL,W,aM) 

^dr  ~  *^dr^  ^2  16) 

6^  =  bJl,L,E,W,aM) 

5^,  =  b,,{l,L,E,W,c,^) 

E  is  Young’s  modulus,  and  W  is  the  specimen  width. 

Note,  that  within  this  sub-subsection,  and  in  every  numerical  procedure  of  this  thesis,  when 
is  used,  it  is  a  negative  quantity.  With  respect  to  this  sub-subsection,  being  negative  will 

result  in  and  6^^  being  negative. 

= -abs(aj^)  (for  all  calculations  in  this  thesis)  (2.17) 

2.2AA  Calculation  of  the  energy  release  rate 

Once  the  CODs  have  been  determined,  the  energy  release  rate  (7j)  is 

calculated  using  equation  2.3.  7j  is  equal  to  the  amount  that  the  system  potential  energy 

decreases  when  the  crack  extends  a  distance  such  that  one  area  unit  of  new  crack  surfaces  is 
generated. 

2.2.4.5  Calculation  of  the  driving  forces 

The  driving  forces  for  crack  extension  X[  and  crack  layer  extension  are 

calculated  using  equations  2.2.  It  is  mathematically  possible  for  either  of  them  to  be  negative, 
which  would  cause  the  crack  length  and/or  the  crack  layer  length  to  decrease.  This  would  be 
physically  impossible,  so  whenever  one  of  the  quantities  is  calculated  and  found  to  be  negative,  it 
is  set  equal  to  zero. 
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2.2.4.6  Calculation  of  the  growth  increments  and  new  values  for  the  crack  and  crack 
layer  lengths 

The  growth  rates  for  the  time  step  are  calculated  using  equation  2.1.  Using 
these  rates,  and  the  time  increment  for  the  time  step  (At),  the  growth  increments  for  the  crack 
( A/ )  and  the  process  zone  ( AL )  are  calculated  by  the  forward  Euler  method. 

A/  =  /-At  (a)  (2.18) 

AL  =  LAt  (b) 

The  new  values  of  the  crack  length  and  crack  layer  length  (L^^^)  are 

^new  “  ^old  1 0'l 

2.2.4.7  Update  of  program  arrays 

The  program  arrays  are  now  be  updated,  and  then  the  calculation  loop  of 
subsection  2.2.4  is  restarted  from  the  beginning. 

l[n]  =  (l^,J,  =  l[n-l]  +  Al 

L[n]  =  (L„,^),  =  L[n  -  1]  +  AL  (2-20) 

t[n]  =  t[n  -  1]  +  At 

2.2.5  Calculation  of  equilibrium  crack  layer  length 

During  each  time  step,  after  the  candidate  new  values  of  the  crack  length  {1^^^^) ,  and 

the  crack  layer  length  (L^^^)  are  calculated  (sub-subsection  2. 2.4.6),  some  conditions  are  tested 

for  using  those  values.  If  one  or  more  of  the  conditions  are  true,  it  may  be  necessary  to  repeat  the 
time  step  using  a  smaller  time  increment.  Three  of  the  conditions  require  the  value  of  the 
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equilibrium  crack  layer  length  .  So  it  is  necessary  to  calculate  .  Always,  the  value  of 
is  the  upper  limit  for  the  value  of  L  (I  <L<  )  (Figure  5). 


Figure  5.  Crack  layer  configuration 


For  a  fixed  value  of  / ,  is  the  value  of  L  which  causes  the  driving  force  for  crack  layer 
extension  to  equal  zero.  From  equation  2.2(b),  =  0  can  be  true  in  two  ways. 


+  =  0  (a) 

K^  +  {l+2y)K^^  =  0  (b)  ■ 


(2.21) 


In  equations  2.21,  (b)  gives  the  shorter  (and  stable)  value  of  L.^  (0<y<oo)  (Figure  6).  (In 

equations  2.21,  (b)’s  value  of  is  shorter  because  always.),  and 

because  increases  from  zero  as  L  increases  from  /.)  Therefore,  while  holding  all  other 
parameters  fixed,  equals  the  value  of  L  which  makes 
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K^  +  (l+2y)Kj^  =  0.  (2.22) 

Therefore,  define 

fiy,  W,  (5jx),l,L)  =  K^  +  {l+  2y)Kj^ .  (2.23) 

Define 

g{L)=  value  of /(L, ...),  for  fixed  values  of  Y,  W,  a^(x), /.  (2.24) 

Then, 

W,  a^(x),  /)  =  ,  L  so  that  g(L)  =  0 .  (2.25) 

For  a  particular  simulation,  y,  W ,  and  are  constants;  and  ctoo(.^)  is  a  prescribed  function. 
Therefore,  for  any  value  of  I  (such  as  the  new  value  of  /  calculated  during  a  time  step),  can 
be  calculated.  The  value  of  L  which  makes  g{L)  =  0,  can  be  found  by  any  method  for 


determining  the  root  of  a  nonlinear  function  (g{L)).  The  root  is  located  within  the  interval 
I  <  root  <  W . 


a)  Showing  maximum  possible  value 


b)  Showing  lower  limit 


Figure  6.  Equilibrium  crack  layer  lengths 
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2.2.6  Enforcement  of  constraints  on  the  crack  and  crack  layer  len2ths 

The  test  conditions  referred  to  in  subsection  2.2.5  are 


'  >  T 

new  new 

(a) 

^new  ^  ^eq(new) 

(b) 

^eq(old)  “  ^old  ^  ^ 

(c) 

^eq(old) 

^new  “  '^eq(new) 

— j — -  <  e 

(  now/ ^ 

(d) 

(2.26) 


e  is  a  prescribed  small  number,  (e  «  1 ). 

If  (a)  and  (b)  are  both  false,  then  the  time  step  continues  (sub-subsection  2.2.4.7). 
If  (a)  is  true,  then  the  time  step  is  repeated^ with  an  adjusted  time  increment. 


Ar  =  (At) 


-L 


^old  ~  ^old 
new)  (•^old 


(Figure  7) 


(2.27) 


If  (b)  and  (c)  are  true,  then  the  value  of  is  adjusted,  and  then  the  time  step  continues  (sub¬ 
subsection  2.2.4.7). 

^new  =  ^e^(new)  (2.28) 

If  (b)  and  (d)  are  true,  and  (c)  is  false,  the  value  of  is  also  adjusted  according  to  equation 
2.28,  and  then  the  time  step  continues  (sub-subsection  2. 2.4.7). 

If  (b)  is  true,  and  (c)  and  (d)  are  false,  then  the  time  step  is  repeated,  with  an  adjusted  time 
increment. 


(l)Within  section  2.2,  whenever  a  time  step  is  repeated,  the  repetition  begins  with  sub-subsection  2.2.4.6. 
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At  =  (At) 


•^eq(old)  ^old 


(^new  -^eqCnew))  (^eq(old)  ^old) 


(2.29) 


During  the  first  attempt  for  all  time  steps,  the  time  increment  (At)  is  set  equal  to  the 
prescribed  initial  time  increment  (A^q)  .  If  the  first  attempt  fails,  then  Ar  will  be  adjusted,  so  that 

during  the  second  attempt  At  <  At^ .  There  is  no  limit  on  the  number  of  attempts  that  a  time  step 


may  require,  but  At  =  Afp  only  on  the  first  attempt. 


Figure  7.  Time  increment  adjustment 


2.2.7  Prevention  of  crack  extension  into  undegraded  material 

Another  program  array  is  introduced  in  this  subsection,  for  the  driving  force  for  crack 
extension,  . 

Define  Ar^^^Q^ ,  so  that 
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2Y(A/p„o,).0.  (2.30) 

A^pz(O)  is  elapsed  time  for  the  resistance  to  crack  extension  of  a  material  subjected  to  the 
drawing  stress,  to  reach  zero  (or  at  least  to  reach  a  value  close  to  zero).  (Equation  2.30  uses 
“  =  instead  of  because  some  material  surface  energy  functions  (2y(...))  decrease 
exponentially.  For  an  exponentially  decreasing  {...),  only  2y{°°)  =  0.)  The  time  increment 
for  the  current  time  step  is  A/ .  If 

— — —  «  1 ,  (in  practice,  is  always  true  (Figure  8.))  (2.3 1) 

then  it  may  be  numerically  possible  for  the  crack  to  extend  into  undegraded  material.  But,  it  is 
normally  physically  impossible  for  the  crack  to  extend  into  undegraded  material.  The  purpose  of 
the  procedure  within  this  subsection  is  to  prevent  the  physically  excluded  crack  extension. 
Procedure:  If  at  the  end  of  the  n  th  time  step  (after  the  calculation  loop  of  subsection  2.2.4  has 
been  successfully  completed). 


1]>0 
X,[n]  =  0 


(if  both  are  true) 


(2.32) 


then,  the  (ordered)  array  /[...]  is  searched  backwards  from  l[n],  to  find  /[/],  so  that 
l[i]  =  /[/  -  1  ] .  What  makes  it  possible  to  find  /[(] ,  so  that  /[(]  =  /[/  -  1  ] ,  is  that  a  crack  may 
sometimes  not  grow  for  a  sequence  of  time  steps.  If  one  or  both  conditions  of  equation  2.32  are 
false,  then  the  n  th  time  step  is  accepted  as  valid.  If  equations  2.32  are  true,  and 

^  >  ( 1  +  e) ,  (e  is  a  prescribed  small  number,  (e  «  1 ).)  (2.33) 

L[i]-l[n-l] 

is  true,  then  the  time  step  is  repeated  with  a  smaller  time  increment. 
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4 


t 


Figure  8.  Penetration  of  crack  into  undegraded  material 


At 


(AO 


L[i]-l[n-l] 

l[n]-l[n-l] 


(1+e). 


Otherwise,  the  n  th  time  step  is  accepted  as  valid. 

2.2.8  Time  jumps  to  decrease  required  computer  time 
Recall  equation  2.1 


/  =  jtj  • 
ij  “  ^2 


(2.34) 


and  equation  2.18 


(a) 

(b) 


(2.35) 
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Al  =  i-At  (a)  ^  (2.36) 

AL  =  LAt  (b) 

These  equations  indicate  that  if  the  calculated  values  of  the  driving  forces  for  crack  length 
extension  (X^)  and  crack  layer  length  extension  (X^)  are  negative  (in  which  case  they  will  be  set 
equal  to  zero),  or  zero,  then  the  crack  and  crack  layer  lengths  cannot  increase.  The  parameters  of 
Xi  and  X^  are 

X;  =  X;(/,  L,  E,  y,  W,  At^^)  (a)  37) 

X^  =  X,(/,  L,  E,  y,  W,  a^(x))  (b)  ' 

For  fixed  values  of  /  and  L ,  the  only  parameter  in  equations  2.37  which  changes  with  time,  is 
Atp^  (the  elapsed  time  since  the  material  at  the  current  location  of  the  crack  tip  (at  I )  was  drawn, 

i.e.,  since  that  material  became  part  of  the  process  zone). 

If  the  current  value  of  X^  is  zero,  and  the  current  value  of  X^  is  very  small  (the  calculated 

value  of  X^  will  always  be  positive,  but  it  can  approach  zero),  then  the  system  will  be  effectively 
frozen  until  enough  time  has  passed,  so  that  At„.  has  increased  and  become  greater  than  some 

pz 

critical  value  (Ar^^^^^^).  Recall  that  Xy  =  J  y-2y{Atp^ .  will  make  the  equality 

/j  =  2y{Atp^(^^^^)  true.  The  maximum  time  increment  used  within  a  simulation  is  A^q  (the  user 

prescribed  time  increment,  which  is  used  during  the  first  attempt  at  executing  every  time  step). 
Depending  on  the  particular  set  of  input  parameters,  it  can  be  possible  that  during  a  simulation 
there  will  be  multiple  instances  when  A^q  «  (and,  X^  =  0  and  X^^ «  1 ).  In  these 

instances,  many  time  steps  will  be  taken  during  which  the  only  system  parameter  which  changes 
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significantly  is  time  (t).  This  can  result  in  a  simulation  which  requires  a  great  amount  of  computer 
time.  The  purpose  of  the  procedure  within  this  subsection  is  to,  when  possible,  instantly  increment 
the  time  by  an  amount  (Atj),  which  causes  the  equality  to  be  true  (Figure  9.). 


Procedure:  At  the  end  of  the  nth  time  step  (if  the  procedure  of  subsection  2.2.7  does  not  cause 
the  time  step  to  be  repeated),  the  following  conditions  are  tested  for 

X,[n]  =  0 

L  [n]-L[n]  •  (e  is  a  prescribed  small  number,  (£«  1 ).)  (2.38) 

_ i ^  r* 
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If  both  of  conditions  2.38  are  true,  then  a  time  jump  can  be  made.  (If  a  time  jump  cannot  be  made, 
then  the  next  time  step  begins  normally.)  Define 

X;(Ati)Hyj-27(Ar,).  (2.39) 

The  root  (At^)  of  function  2.39,  is  found  by  any  method  for  finding  the  root  of  a  nonlinear 
equation.  The  lower  limit  of  the  interval  which  contains  Ar^  is  the  current  value  of  .  The 
upper  limit  is  Ar^^^Qj  (see  equation  2.30). 

<  ^^pz(0)  (2-40) 

Define 

Atj^At^-Atp^.  (2.41) 

Now  the  time  step  counter  is  incremented  and  then  the  array  values  are  set  for  this  new  time  step, 

n  =  «  +  1 

l[n]  =  l[n-\]  (2.42) 

L[n]  =  L[n-  1] 

t[n]  =  t[n  -  1]  +  Atj 

The  calculation  loop  of  subsection  2.2.4  is  now  restarted  from  the  beginning. 

2.2.9  Determination  of  specimen  failure 

The  criterion  for  specimen  failure,  is  that  L ,  the  crack  layer  length,  equals  (approaches 
very  close  to)  W,  the  specimen  width.  Equivalently,  the  specimen  fails  when  all  (almost  all)  of  the 
material  ahead  of  the  crack  tip  (at  / )  is  drawn  (when  the  process  zone  reaches  (almost  reaches)  the 


far  edge  of  the  specimen). 
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Procedure:  During  each  time  step,  after  the  candidate  new  values  of  the  crack  length  ,  and 
the  crack  layer  length  are  calculated  (sub-subsection  2.2.4.6),  and  immediately  before  the 

procedure  of  subsection  2.2.6  is  performed,  the  following  condition  is  tested  for. 

>  W  (2.43) 

If  condition  2.43  is  false,  then  the  time  step  continues  normally  with  subsection  2.2.6.  If  condition 
2.43  is  true,  then  the  time  increment  is  adjusted,  and  the  time  step  repeated. 

.  At 

Ar  =  -  (2.44) 

If  during  the  nth  time  step,  condition  2.43  becomes  true  for  the  sixth  time,  the  nth  time  step  is 
abandoned,  the  specimen  is  presumed  to  have  failed,  and  the  simulation  is  ended.  The  lifetime  (or 
time  to  failure)  (ty )  of  the  specimen  is 

ty=t[n-l].  (2.45) 

Using  this  method  to  determine  specimen  failure,  insures  that  the  final  value  of  the  crack  layer 
length  (Ly)  satisfies 

W-Lf 

<<  1  •  (2.46) 

2.3  Algorithm  for  SEN  specimen  sub.!ected  to  load  control 

This  section  will  specify  the  changes  which  are  made  to  the  algorithm  of  section  2.2 
(constant  remote  load)  when  the  remote  load  varies  with  time,  i.e.,  when  the  remote  load  function 


changes  from  cy^(x)  to  a„(x,  t). 
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A  new  sub-subsection  (2.2.4.1(a))  is  added,  in  which  the  remote  load  function  (ct^(x,  t))  is 
updated  at  the  beginning  of  each  time  step  for  the  new  time.  For  instance,  at  the  start  of  the  first 
time  step,  the  function  would  be  o^(j:,  0) . 

The  energy  release  rate  (J^ )  is  now  (explicitly)  time  dependent.  It  seems  logical  that  the  time 

jumps  of  subsection  2.2.8  could  be  used  in  this  algorithm  (load  control)  by  modifying  equation 
2.39  to  become 

)  =  7i (t)  -  2y ( At  1 ) .  (t  =  tp^  + At,)  (2.47) 

Doing  this  would  be  wrong.  The  algorithm  for  load  control  cannot  allow  time  jumps.  For  a  time 
varying  remote  load  function,  equations  2.37  become 

Xi  =  X;(/,  L,  E,  y,  W,  <5^{x,  t),  At^^)  (a)  4g) 

X^  =  X,il,L,E,y,c^^,W,oJx,t))  (b) 

The  driving  force  for  crack  layer  extension  (X^ )  is  now  time  dependent.  This  means  that  if  a  time 

jump  was  made,  the  crack  layer  length  (L)  would  increase.  It  would  be  necessary  to  be  able  to 
calculate  the  size  of  the  increase  in  L,  in  order  to  be  able  to  calculate  J,  in  equation  2.47. 

However,  calculating  the  value  that  L  would  have  after  a  time  jump  was  made  into  the  future  is 
impossible.  The  only  way  to  calculate  the  value  of  L  for  a  time  t ,  is  to  perform  the  calculation 
loop  of  subsection  2.2.4  until  time  t  is  reached;  i.e.,  there  is  no  equation  which  will  provide  a 
future  value  of  L .  Therefore,  every  advance  in  time  must  be  made  using  subsection  2.2.4,  and  the 
time  jump  procedure  of  subsection  2.2.8  cannot  in  any  way  be  used  in  this  (load  control) 
algorithm. 
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2.4  Assumption  of  simplified  crack  layer  model  for  discontinuous  growth 


The  material  surface  energy  function  {2y{Atp^))  usually  has  one  of  two  forms. 


2Y(Af  )  =  2Yoe 


(exponential)  (2.49) 


or 


2Yo(l-fco^tp,)<0: 


2Y(At^,)  =  2Yo(l-/coAV) 

2y{Atp^)  =  0 


(linear) 


(2.50) 


2Yo  ,  a  material  constant,  is  the  surface  energy  for  drawn  material  with  no  degradation.  At^^  is  the 

elapsed  time  since  the  material  at  the  current  location  of  the  crack  tip  (at  I )  was  drawn.  kQ ,  a 

material  constant,  is  the  degradation  coefficient.  (The  information  in  this  (first)  paragraph  will  be 
used  later  in  this  section.) 

In  industrial  applications,  plastics  are  usually  subjected  to  a  remote  load  function  (a„(A:)), 
whose  average  value  (<7^^gygP  is  small  compared  to  the  material  drawing  stress 
Laboratory  test  specimens  subjected  to  a  loading  in  which  is  small,  have  been  observed  to 

exhibit  discontinuous  crack  and  crack  layer  growth.  By  “discontinuous  growth”  is  meant;  a  long 
time  interval  during  which  the  crack  is  stationary  and  the  crack  layer  is  almost  stationary,  is 
followed  by  a  very  short  time  interval  during  which  the  crack  and  crack  layer  grow.  This  cycle 
repeats  until  the  specimen  has  almost  failed,  when  the  growth  becomes  continuous.  During  the 
discontinuous  part  of  the  growth,  each  new  stationary  location  of  the  crack  tip  (at  /)  is  near  the 
previous  (almost)  stationary  location  of  the  front  of  the  crack  layer  (at  L).  This  type  of  crack  and 
crack  layer  growth  is  referred  to  as  “discontinuous”,  “stepwise”  or  “staircase”  growth. 
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Discontinuous  growth  similar  to  what  has  been  experimentally  observed,  can  be  numerically 
simulated  if  the  set  of  input  parameters  is  adjusted  so  that  two  conditions  are  satisfied.  Condition 
1  is  that  the  value  of  the  energy  release  rate  at  the  start  of  the  simulation  )’  than  the 

material’s  initial  surface  energy  (2Yo)-  Whenever  J ^  <  2y{Atp^) ,  growth  is  discontinuous.  This  is 
because  the  driving  force  for  crack  extension,  Xf ,  equals  J ^-2y {At p^) .  When  is  calculated  to 
be  negative  (in  which  case  it  is  set  equal  to  zero),  or  zero,  then  time  passes  during  which  the  crack 
does  not  grow.  During  a  creep  (constant  load)  simulation  (or  test),  7|  is  an  increasing  function,  so 

that  if  and  when  7|  becomes  greater  than  2yQ  (the  maximum  possible  value  of  2y{Atp^)),  the 

remainder  of  the  test  will  exhibit  continuous  growth,  i.e.,  the  crack  will  then  always  be  growing. 
Satisfaction  of  condition  one  guarantees  that  the  growth  is  discontinuous  (at  least  initially). 

The  second  condition  required  to  simulate  experimentally  observed  discontinuous  growth,  is 

that 

k^,k2»kQ.  (2.51) 

itj ,  a  material  constant,  is  the  crack  length  kinetic  coefficient.  k2 ,  a  material  constant,  is  the  crack 
layer  length  kinetic  coefficient,  k^  is  the  degradation  coefficient  mentioned  in  this  section’s  first 
paragraph.  Recall  equation  2.1, 


l  =  k^Xi 

(a) 

L  =  ^2  • 

(b) 

(2.52) 


Satisfaction  of  restriction  2.5 1  (condition  2)  guarantees  that  the  time  intervals  between  growth  are 
much  larger  than  the  time  intervals  of  growth. 
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During  a  numerical  simulation  of  discontinuous  growth,  a  cycle  is  repeated  in  which:  (a)  the 
crack  and  crack  layer  simultaneously  begin  to  grow  very  rapidly;  (b)  the  crack  stops  growing 
while  the  crack  layer  continues  to  grow,  but  very  slowly  and  for  a  large  time  interval.  Each  time 
the  crack  stops  growing,  the  crack  layer  will  have  almost  reached  the  equilibrium  crack  layer 
length  The  crack  layer  length  (L)  will  then  asymptotically  approach  the  equilibrium  crack 

layer  length,  until  another  period  of  rapid  growth  begins.  (L  approaches  asymptotically, 

c  y 

because  the  driving  force  for  crack  layer  length  extension  X^(L, ...)  is  a  decreasing  function  of 
L ,  having  maximum  value  when  L  =  I,  and  minimum  value  (zero)  when  L  =  L . _ .)  A  natural 

simplification  for  modeling  discontinuous  growth,  is  to  assume  that  the  value  of  L  is  always  equal 
to  the  value  of  L  . 

L(t, ...)  =  Lg^(t,  ...)  (simplification  for  modeling  discontinuous  growth)  (2.53) 

When  this  simplification  is  made,  the  parameters  (for  a  creep  specimen)  of  the  crack  layer  length 
(L)  become 

L  =  L^qiy,  w,  a^(x),  /) .  (2.54) 

The  only  parameter  in  equation  2.54  which  is  not  constant  during  a  simulation  is  the  crack  length 
(/).  Therefore,  for  a  particular  simulation,  the  crack  layer  length  is  a  function  of  only  the  crack 
length  (L  =  Lg^(l) ).  So,  the  crack  layer  grows  when  and  only  when  the  crack  grows.  Also,  all  of 

the  process  zone  material  (in  front  of  the  crack  tip)  will  be  equally  degraded,  because  the  crack 
layer  cannot  grow  while  the  crack  is  stationary.  This  means  that  each  instance  of  crack  growth 
consists  of  an  instantaneous  jump  to  the  end  (at  of  the  process  zone.  The  growth  cycle 
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consists  of:  (a)  time  passing  until  the  driving  force  for  crack  extension  (Xj ,  =  /j  -  2y(Atp^) )  is 
positive  (until  >  ly{Atp^) );  (b)  the  crack  tip  (/^;^)  instantaneously  jumping  to  the  end  of  the 
process  zone  while  the  front  of  the  crack  layer  instantaneously  jumps  to  the 

new  position  ■ 

Using  the  simplification  of  equation  2.53  reduces  the  complexity  of  the  numerical  algorithms 
for  modeling  discontinuous  growth,  and  results  in  simulations  which  require  relatively  little 
computer  time.  The  lifetime  (or  time  to  failure)  (/y)  of  a  specimen  becomes  slightly  (but  not 
significantly)  shorter.  So  a  simulation  based  on  the  simplified  CLM  gives  a  somewhat  more 
conservative  value  for  tj-  than  the  same  simulation  based  on  the  (complete)  CLM  does  (Figure 
10.). 


Lj  j/j 

hA 


Figure  10.  Simplified  model  of  discontinuous  crack  growth 
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2.5  Simplified  algorithm  for  SEN  specimen  subjected  to  creep  loading 

This  algorithm  uses  the  simplifying  assumption  of  section  2.4,  i.e.,  that  the  crack  layer 
length  (L),  is  always  equal  to  the  equilibrium  crack  layer  length  {L^^). 

2.5.1  Determination  of  remote  load  function 
See  subsection  2.2.2. 

2.5.2  Program  Arrays 

The  initial  conditions  (of  subsection  2.2.1)  now  become 
/[O]  =  0 

l[0]  =  Iq  .  (2.55) 

L[0]  =  L,^(/o) 

Recall  from  equation  2.25  that 

^eq  =  /)  .  (2.56) 

With  I  =  Iq,  all  parameters  of  equation  2.56  are  known,  and  L[0]  can  be  calculated  according  to 
the  procedure  in  subsection  2.2.5. 

2.5.3  Calculation  loop  for  time  evolution  of  crack  and  crack  layer 

During  odd  time  steps  (1,3,  ...),  time  (t)  passes,  while  the  crack  length  (0  and  crack 
layer  length  (L)  are  stationary.  During  even  time  steps  (2, 4, ...),  the  crack  length  (/ )  and  the  crack 

layer  length  (L)  both  grow,  while  the  time  is  stationary.  (“Time  step”  is  actually  a  misnomer  for 
the  even  “time  steps”,  because  for  them  the  time  increment  is  zero.) 

2.5.3. 1  Increment  of  time  step  counter 


See  sub-subsection  2.2.4. 1. 
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2.5.3.2  For  odd  time  steps 

2.5.3.2.1  Calculation  of  CODs 

The  crack  opening  displacements  (CODs),  5^,  and  6^^  are  calculated 

according  to  sub-subsection  2.2.4.3. 

2.5.3.2.2  Calculation  of  the  energy  release  rate 
See  sub- subsection  2.2.4.4. 

2.5.3.2.3  Calculation  of  time  increment 

The  procedure  in  this  sub-sub-subsection  is  similar  to  the  procedure  in 

subsection  2.2.8. 

Define 

/(Ar)  =  /i-2Y(Ar).  (2.57) 

The  time  increment  is  the  root  (Ar^)  of  function  2.57,  which  can  be  found  by  any  method  for 
finding  the  root  of  a  nonlinear  function. 

0  <  Ar^  <  Atp^^Qj .  (2.58) 

Afp^^Q^  is  defined  near  equation  2.30.  According  to  the  assumption  (L  =  of  section  2.4,  the 

crack  tip  (/)  and  the  front  of  the  crack  layer  (L)  have  just  jumped  to  new  positions  (except  for 
during  the  first  time  step,  when  t  =  0),  and  all  of  the  process  zone  material  ahead  of  the  crack  tip 
is  now  beginning  to  degrade.  That  is  why  the  lower  limit  of  interval  2.58  is  zero. 

As  the  simulation  proceeds,  and  the  crack  and  crack  layer  lengths  grow,  the  value  of  the 
energy  release  rate  (J^ )  increases,  and  the  elapsed  (degradation)  time  (Ar^)  until  the  crack  and 


crack  layer  jump  to  new  positions,  decreases. 
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2.5.3.2.4  Update  ofprosram  array 

The  time  array  is  now  updated,  and  then  the  calculation  loop  of  subsection 
2.5.3  is  restarted  from  the  beginning. 


t[n]  =  t[n-\]  +  At^  (2.59) 

2.5.3.3  For  even  time  steps 

2.5.3.3.1  Calculation  of  the  new  values  for  the  crack  and  crack  layer  lengths 

The  crack  tip  jumps  to  the  end  of  the  current  process  zone,  and  the  front  of 
the  crack  layer  jumps  to  the  equilibrium  length  for  the  new  position  of  the  crack  tip. 


^new  ^old 
^new  ~  ^eq(^new) 


(2.60) 


^eqi^new^  calculated  according  to  the  procedure  in  subsection  2.2.5. 

2.5.3.3.2  Update  of  program  arrays 

The  crack  and  crack  layer  length  arrays  are  now  updated,  and  then  the 
calculation  loop  of  subsection  2.5.3  is  restarted  from  the  beginning. 


=  /„ew 

=  4ew 


(2.61) 


2.5.4  Determination  of  specimen  failure 

When  becomes  equal  to  IJq  (the  maximum  value  for  the  growth 

changes  from  discontinuous  to  continuous  (i.e.,  for  the  remainder  of  the  simulation  no  time  will 
pass  without  crack  and  crack  layer  growth).  For  the  simplified  CLM,  increments  in  crack  and 
crack  layer  growth  are  accomplished  in  zero  time.  Since,  after  the  onset  of  continuous  growth. 
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growth  must  always  be  occurring,  and  since  all  growth  must  happen  instantaneously,  once 
continuous  growth  begins  time  cannot  increase.  Therefore,  the  onset  of  continuous  growth  must 
be  the  indicator  of  specimen  failure.  So,  the  criterion  for  specimen  failure  is  =  2Yo  (Figure 
11). 

Procedure;  After  each  odd  time  step  is  made,  the  condition  7^  >  2Yo  is  tested  for.  If  the  condition 
is  true,  then  the  simulation  is  ended,  and  the  lifetime  (t^)  (or  time  to  failure)  for  the  specimen  is 
tf  =  tin] .  (2.62) 


Figure  1 1 .  Determination  of  lifetime 


2.6  Simplified  algorithm  for  SEN  specimen  sub.tected  to  load  control 

This  section  will  specify  the  changes  which  are  made  to  the  algorithm  of  section  2.5 
(constant  remote  load)  when  the  remote  load  varies  with  time,  i.e.,  when  the  remote  load  function 
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changes  from  a^{x)  to  a^(;r,  t) .  The  algorithm  uses  the  simplifying  assumption  of  section  2.4, 


i.e.,  that  the  crack  layer  length  (L),  is  always  equal  to  the  equilibrium  crack  layer  length  (L^^). 

A  new  sub-subsection  (2.5.3.1(a))  is  added,  in  which  the  remote  load  function  (g<„(;c,  t))  is 

updated  at  the  beginning  of  each  new  time  step,  whenever  the  time  changes. 

2.6.1  Program  Arrays 

This  subsection  (2.6.1)  replaces  subsection  2.5.2. 

The  initial  conditions  are  the  same  as  specified  in  subsection  2.5.2.  The  parameters  of  L 

eq 

(see  equation  2.25)  become 

^eq  =  Legiy,  1^,  oM,  0,  0  •  (2.63) 


With  I  —  Iq  and  r  —  0,  all  parameters  of  equation  2.63  are  known,  and  L[0]  can  be  calculated 


according  to  the  procedure  in  subsection  2.2.5. 

2.6.2  Addendum  for  calculation  of  time  increment 

This  subsection  (2.6.2)  modifies  sub-sub-subsection  2.5. 3.2.3. 

Define 


^/(^i)  =  ^i(^i)-2y(Ar^^).  (2.64) 

Equation  2.64  replaces  equation  2.57. 

The  root  (t^)  of  function  2.64,  is  found  by  any  method  for  finding  the  root  of  a  nonlinear  equation. 
The  root  search  interval  is 


t<t^<t  + 


^^pz(0)  • 


(2.65) 


Equation  2.65  replaces  equation  2.58.  Each  time  function  2.64  is  evaluated,  the  surface  energy 
C^yi^t pz))  is  calculated.  The  remote  load  function  (C7^(.x,  tj))  is  updated.  The  crack  layer 
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length  (L  =  a^(x,  tj), ...))  is  calculated.  The  CODs,  5^  =  5^(/,  L,  E,  W,  <5^{x,  tj))  and 

6^^  =  L,  W,  are  calculated.  The  energy  release  rate 

(7j  =  (5^  +  ( 1  +  y)S^r)(“^Jr))  (equation  2.3)  is  calculated. 

Recall  that  in  section  2.3  (complete  algorithm  for  load  control)  the  size  of  a  time  jump  could 
not  be  calculated,  because  the  value  of  the  crack  layer  length  (L)  was  history  dependent.  Within 
this  algorithm  (section  2.6,  simplified  algorithm  for  load  control),  the  size  of  a  time  jump  (time 
increment)  can  be  calculated,  because  the  simplifying  assumption  L  =  L  is  used.  The  value  of 

the  equilibrium  crack  layer  length  (L  )  is  not  history  dependent  (Figure  12.). 


Figure  12.  Time  jump  under  increasing  load 
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2.7  Basis  for  simplified  algorithm  for  SEN  specimen  under  displacement 

CONTROL 

When  the  load  is  a  prescribed  function  of  time  (as  is  true  for  specimens  subjected  to  creep  or 
load  control),  the  only  relevant  parameter  of  specimen  geometry  is  the  specimen  width  (W)  (the 
specimen  thickness  (5)  is  not  relevant,  (subsection  2.2.2)).  When  displacement  is  a  prescribed 
function  of  time  (as  is  true  for  specimens  subjected  to  displacement  control),  an  additional 
parameter  of  specimen  geometry  is  required,  the  specimen  height  (2H).  2H  is  the  initial 
(unstrained)  distance  between  the  two  grips  which  hold  the  specimen.  H  is  the  initial  distance 
from  the  centerline  of  the  crack  to  either  of  the  two  grips. 

For  the  displacement  control  mode  of  loading,  the  distance  {2v{x,  /))  between  the  two  grips 
is  prescribed  as  a  function  of  time.  The  remote  load  function  (a^(A:,  t) ,  the  traction  at  either  grip) 

must  be  calculated.  An  approximation  that  will  be  used  in  order  to  make  the  calculations  possible, 
is  that 

0  =  •  (approximation  used  for  displacement  control  algorithm)(2.66) 

In  other  words,  at  any  time  t ,  the  remote  load  stress  will  be  assumed  to  be  uniform.  In  order  to 
calculate  o„(0 ,  the  displacement  which  is  conjugate  to  it  must  be  calculated.  That  displacement 
is 

a:=  W 

1  (2.67) 

x^O 

V^(t)  is  the  average  value  of  the  prescribed  displacement  {v{x,  t))  at  the  top  of  the  specimen,  at 
time  t  (the  x-axis  is  located  at  the  centerline  of  the  crack).  (The  average  value  of  the  prescribed 
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displacement  at  the  bottom  of  the  specimen  at  time  t  would  be  —Vp{t) .)  v{x,  t)  can  be  divided 
into  two  parts.  The  first  part  is  due  to,  a^(f)  applied  to  the  solid  (uncracked)  body,  and  is  equal  to 


n(o  = 


E 


(2.68) 


{E  is  Young’s  modulus.)  The  second  part  of  v{x,  t) ,  which  will  be  called  V^{x,  t) ,  is  due  to  the 
tractions  a^(0  and  acting  on  the  top  crack  face  (0  <  o^(t)  <L,  I  <  ^L).  The  material 

of  the  process  zone  (I  <  x  <  L)  is  imagined  to  be  removed,  so  that  the  crack  extends  from  x  =  0 
to  a:  =  L.  The  process  zone  material  is  then  supposed  to  be  replaced  by  tractions  which  act  to 
compress  the  crack,  and  are  equal  in  magnitude  to  the  process  zone  stress  Also,  when  it  is 
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pretended  that  the  process  zone  material  is  removed,  the  specimen  becomes  linearly  elastic. 
According  to  Betti’s  reciprocal  theorem  (Figure  13), 


Figure  13.  Calculation  of  edge  displacement 


^  ^  G  (0^  ~  ^  G  ^  ~  ^ 

1-  J  v^{x,t)dx  =  I  b^{x,G^=\,...)dx  +  -^  J  ?>^{x,o^=\,...)dx.  (2.69) 

x  =  Q  a:  =  0  X  =  l 

(Soo(^>  <^00  =  •••)  is  ths  crack  opening  displacement  at  x-coordinate  x,  due  to  the  remote  load 

function  (stress)  =  1 .)  V’p(r)  can  now  be  written 
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x^W 


J  {V,W  +  V^(x.t))dx. 


(2.70) 


x  =  0 


Therefore, 


x  =  L 


Vp(0  =  ^^^  +  ^^  J  5^(x,a^  =  1,  ...)6?^+ J  5^(x,a^  =  1,  ...Mj^(2.71) 

jc  =  0  X  -  I 


So, 


x-L 

''p<»)-^  J  8.(x,ct,  =  I,  ...)& 

0,(0 - - .  (2.72) 

I +  ^  J  5jx,a^= 

x  =  0 

The  simplified  algorithm  for  displacement  control,  is  like  the  simplified  algorithm  for  load 
control.  At  the  beginning  of  each  time  step,  time  is  incremented.  Then  the  prescribed  remote 
displacement  function  (v(x,t))  is  updated.  From  v{x,t),  V^it)  is  calculated.  Then  using  the 

procedure  of  subsection  2.2.5,  the  crack  layer  length  L  =  L  (simplified  assumption,  section 
2.6)  is  calculated.  Then  using  equation  2.72,  a^(r)  can  be  calculated.  After  is  calculated, 

the  algorithm  proceeds  like  the  simplified  load  control  algorithm.  But  there  is  a  problem  regarding 
the  value  of  L .  The  parameters  of  L  are 

L  =  (2.73) 

From  equation  2.72,  the  parameters  of  c^{t)  are. 


Cf=o(0  =  L,  E,  W,  H,  vix,  /)) . 


(2.74) 
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So,  L  =  L(o^(0,  •■•)  >  (y^{t)  =  <5„{L,  ...),  i.e.,  each  is  a  function  of  the  other.  Therefore, 

L  and  a^(t)  must  be  computed  simultaneously  in  a  cycle  during  which  the  value  of  the  first  is 
used  to  calculate  the  second,  and  then  the  value  of  the  second  is  used  to  calculate  the  first.  The 
cycle  continues  until  the  change  in  L  and  G^(r)  between  successive  cycles  is  very  small. 

The  computation  of  L  and  (T^(0  is  the  main  difficulty  associated  with  the  implementation  of 
the  simplified  displacement  control  algorithm.  It  is  also  the  reason  that  the  complete  algorithm 
{L<  Lg^)  for  displacement  control  does  not  appear  in  this  thesis.  A  simulation  using  the  complete 

algorithm  would  use  a  relatively  great  amount  of  computer  time.  Also,  the  primary  interest  for 
industry  is  discontinuous  crack  and  crack  layer  growth,  which  the  simplified  algorithm  can  model 
adequately. 

2.8  Simplified  algorithm  for  SEN  specimen  sub.tected  to  displacement 

CONTROL 

This  algorithm  uses  the  simplifying  assumption  of  section  2.4,  i.e.,  that  the  crack  layer 
length  (L),  is  always  equal  to  the  equilibrium  crack  layer  length  {L  ). 

2.8.1  Program  arrays 

The  initial  conditions  are  the  same  as  specified  in  subsection  2.5.2,  with  one  addition  specified 
at  the  end  of  this  subsection  (2.8.1).  The  parameters  of  (see  equation  2.25)  become 

L  =  Lg^(Y,  W,  a^(r),  /) .  (2.75) 

I  -  Iq  and  t  =  0.  If  the  value  of  the  prescribed  displacement  function  v{x,  t  =  0)  =  0,  then 
^oo(O)  =  0 ,  all  parameters  of  equation  2.75  are  known,  and  L[0]  can  be  calculated  according  to 
the  procedure  in  subsection  2.2.5.  If  v(x,  t  =  0)9^0,  then  a^(0)  and  L[0]  can  be  calculated 
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according  to  the  procedure  in  sub-subsection  2. 8.2.4.  Since  for  displacement  control  the  remote 
stress  (cr<„(t))  must  be  calculated,  its  values  are  recorded  in  an  array  (a^[i]).  The  additional 

initial  condition  is 

cJO]  =  a^(0).  (2.76) 

2.8.2  Calculation  loop  for  time  evolution  of  crack  and  crack  layer 

2.8.2.1  Increment  of  time  step  counter 
See  sub-subsection  2.2.4. 1. 

2.8.2.2  Increment  of  time 

The  initial  time  step  increment  (A^q,  see  subsection  2.2.3)  is  used  for  every 
time  step.  The  time  {t)  used  during  a  time  step,  is  equal  to  the  previous  time  (t^j^j)  plus  Afp . 

^  =  ^oid  +  ^^o  (2-77) 

2.8.2.3  Update  of  displacement  function 

The  prescribed  displacement  function  (v(x,  t))  is  updated  for  the  current  value  of  the  time 

(0. 

2.8.2.4  Simultaneous  calculation  of  the  remote  load  stress  and  crack  layer  length 
Note  that  in  this  sub-subsection,  <3^{t)  will  be  called  and  Vp(t)  will  be  called 

The  value  for  the  crack  layer  length  (L  =  L^^)  must  satisfy 


W,  a„)  +  ( 1  +  2y)  •  L,  W,  =  0 .  (see  equation  2.22) 


(2.78) 
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In  order  to  insure  convergence  when  calculating  and  L ,  it  is  more  convenient  to  put  equation 
2.72  in  the  form. 


f  x-L 

^  +  ^  J  8jx,/,L,£,W,a^=  IMx  J  =  0.(2.79) 


,  £  21V 

;c  =  0 


,  2W 

J  x  =  l 


The  value  for  the  remote  stress  (a^ )  must  satisfy  equation  2.79.  It  is  more  useful  to  write  that  the 
pair  (a^,  L)  must  satisfy  equations  2.78  and  2.79.  Therefore,  it  is  desired  to  solve  the  system  of 


two  nonlinear  equations  2.78  and  2.79,  for  the  unknown  variables  and  L  (Figure  14). 


Figure  14.  Remote  stress  -  crack  layer  length  relationship 
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Within  the  procedure  to  calculate  the  remote  stress  (o^)  and  the  crack  layer  length  (L),  some 

intermediate  values  will  be  used.  The  intermediate  values  of  will  be  called  , 

,  and  .  The  intermediate  values  of  L  will  be  called  and  . 

It  is  necessary  to  define  some  functions  which  will  be  used  within  the  procedure  to  calculate 
and  L . 


/(oj  = 

function  which  returns  value  of  L,  by  solving  equation  2.78  using 

prescribed  value 

f(L)^ 

function  which  returns  value  of  G^,  by  solving  equation  2.78  using 

prescribed  value  L 

g(<yj  = 

function  which  returns  value  of  L,  by  solving  equation  2.79  using 

prescribed  value  G^ 

g(L)^ 

function  which  returns  value  of  G^,  by  solving  equation  2.79  using 

prescribed  value  L 

mag=  (g^,L) 

function  which  returns  the  vector  magnitude  of  the  pair  (a^,  L) 

The  usual  procedure  to  simultaneously  solve  equations  2.78  and  2.79,  is  to  repeat  a  cycle  in 
which  equation  2.78  is  solved  for  L  using  the  value  of  previously  calculated  by  solving 

equation  2.79,  and  then  equation  2.79  is  solved  for  g^  using  the  value  of  L  previously  calculated 

by  solving  equation  2.78.  But  for  large  values  of  the  prescribed  displacement  rate  (v(x,  t)) ,  this 
method  diverges.  If  is  represented  by  the  vertical  axis,  and  L  is  represented  by  the  horizontal 
axis,  then  when  equations  2.78  and  2.79  are  plotted,  they  appear  in  the  first  quadrant  as  two 
intersecting  curves.  It  is  desired  to  find  the  intersection  point  (L,  G^) .  If  the  sequence  of  points 
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calculated  using  the  procedure  just  described  are  connected,  they  form  a  rectangular  spiral  (Figure 
15(a))  around  the  intersection  point.  For  convergence,  the  spiral  grows  inward,  and  for 
divergence,  it  grows  outward.  The  procedure  which  is  described  below  always  converges.  It  works 
by  always  remaining  on  the  two  upper  branches  of  the  intersecting  curves,  and  monotonically 
moving  towards  the  intersection  point  (Figure  15(b)). 


Figure  15.  Procedures  of  solution  for  two  equation  system 


Procedure: 

First  Fp  is  calculated  using  equation  2.67.  Then  a  numerical  procedure  is  performed  which  is 
presented  next. 

old  r 

<^00  = 


~  1  ]  equals  the  value  of  for  the  previous  time  step) 


63 


1°^^*  =  L[n-1] 


(L[n  -  1]  equals  the  value  of  L  for  the  previous  time  step) 


Oc  =  ) 


L  =  f{c^  ) 


y  old  j  A 
r  new  L  •¥  L 


BEGIN  WOP 


old  .  old  yOld. 

m  =  mag(a^  ,L  ) 


B  ,  r  new . 
<Joo  =  ) 


IF  ^  <^oo^  then 


=  a. 


L"  =  fioj 


j  new  j  A 
r  new  L  +  L 


ELSE 


a„  =  o„ 


,new  ,C 
^new  ^  L  -¥  L 


2 
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END  IF 

L  =  L""'" 

m""'"  =  mag(a^,  L) 

,  .  new  old. 

IF  )  <  e  (e  is  a  prescribed  small  number,  (e  «  1 ).)  THEN  EXIT  LOOP 


END  LOOP 

After  the  loop  is  exited,  the  new  values  of  the  remote  stress  (ct^)  and  crack  layer  length  (L) 
have  been  calculated  and  set. 

2.5.2.5  Determination  of  the  surface  ener^x 

The  surface  energy  is  2Y(Af^^)  (see  sub-subsection  2.2.4.2).  As  indicated  in 
equation  2.12,  At^^  =  For  this  algorithm,  is  the  time  of  the  last  crack  jump,  and  (of 

course)  t  is  the  current  time.  For  the  first  time  step,  =  0. 

2.8.2.6  Calculation  ofCODs 

The  crack  opening  displacements  (CODs),  5^  and  5^^  are  calculated 

according  to  sub-subsection  2.2.4.3  (actually,  according  to  the  appendix). 

2.8.2.7  Calculation  of  the  enerev  release  rate 


See  sub-subsection  2.2.4.4. 
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2.8.2.8  Calculation  of  driving  force  for  crack  length  extension 

The  driving  force  for  crack  length  extension  (X^)  is  calculated  as  specified  in 

sub-subsection  2.2.4.5.  What  happens  next  depends  on  whether  is  zero  or  positive. 

2.8.2.9  If  driving  force  for  crack  leneth  extension  is  zero 

2.8.2.9.1  Update  of  program  arrays 

The  crack  and  crack  layer  lengths,  and  time  arrays  are  now  updated,  and 
then  the  calculation  loop  of  subsection  2.8.2  is  restarted  from  the  beginning. 


l[n]  =  l[n-l] 

L[n]  =  L 
=  <7^ 

r[n]  =  t[n-  1]  +  A/q 


(2.80) 


2.8.2.10  If  driving  force  for  crack  length  extension  is  positive 

If  the  driving  force  for  crack  length  extension  (X^)  is  positive,  then  the  crack 

tip  (/)  will  jump  ahead  to  a  new  position.  It  is  necessary  to  calculate  the  exact  time  (tj)  when 

X[  =  0,  and  the  crack  jumps. 

2.8.2.10.1  Calculation  of  crack  length  jump  time 
Define 

=  h-tpf).  (2.81) 

The  root  {tj)  of  function  2.81,  is  found  by  any  method  for  finding  the  root  of  a  nonlinear  equation. 
The  root  search  interval  is 


^old  <f<t- 


(see  equation  2.77) 


(2.82) 
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In  order  to  determine  the  value  of  function  2.81  for  a  particular  value  ,  t  (sub-subsection 
2.8.2. 2)  is  set  equal  to  .  Then  the  procedures  in  sub-subsections  2.8.2.3-2.8.2.8  are  performed. 
After  the  time  {t^)  when  the  crack  tip  jumps  is  found,  the  time  see  sub-subsection  2. 8.2.5) 
when  degradation  of  the  process  zone  begins  is  reset, 

tpz  =  0-  (2-83) 

Now,  because  the  crack  length  is  making  a  jump,  an  extra  time  step  must  be  taken  before  the 
calculation  loop  of  subsection  2.8.2  begins  again.  During  the  current  time  step,  time  jumps  from 
told  ( t[n -  1  ] )  to  tj,  the  crack  layer  length  grows  from  L[n-l]  to  the  new  calculated  value  L , 

the  remote  stress  grows  from  a^[n  -  1  ]  to  the  new  calculated  value  ,  and  the  crack  length  (/) 
remains  stationary.  During  the  extra  time  step  (starting  with  sub-sub-subsection  2.8.2.10.3)  the 
time  step  counter  is  incremented  (n  =  n  +  1 ),  and  then  the  crack  length  jumps  from  /[n  -  1]  to 
L[n  -  1  ] ,  while  L  and  jump  to  new  values  which  are  calculated  using  the  new  value  of  the 

crack  length.  During  the  extra  “time  step”,  the  time  (t)  and  therefore  the  prescribed  displacement 
(v(x,  r))  remain  stationary.  For  the  extra  time  step,  the  jump  in  the  remote  stress  (a^ )  is  negative, 

because  the  displacement  remains  fixed  while  the  specimen  stiffness  decreases  (because  the  crack 
and  crack  layer  lengths  have  increased). 

2.8.2.10.2  Update  of  program  arrays  for  time  rump 

The  crack  length,  crack  layer  length,  remote  stress  and  time  arrays  are 


updated. 
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l[n]  =  l[n-l] 

~  ^  (L  and  a  are  the  values  calculated  in  sub-subsection  2.8.2.4.)(2.84) 

cr«.[«]  = 

t[n]  =  tj 

2.8.2.10.3  Increment  o  f  time  step  counter  for  extra  time  step 
See  sub-subsection  2.2.4. 1. 

2.8.2.10.4  Simultaneous  calculation  of  the  remote  load  stress  and  crack  layer  length 
for  the  extra  time  step 

The  crack  length  (/)  is  set  equal  to  the  current  value  of  the  crack  layer 
length  (L[n  -  1  ] ),  and  then  new  values  of  the  remote  stress  (a^ )  and  crack  layer  length  (L)  are 


calculated  using  the  procedure  of  sub-subsection  2.8.2.4. 

2.8.2.10.5  Update  of  program  arrays  for  crack  length  lump 

The  crack  length,  crack  layer  length,  remote  stress  and  time  arrays  are 
updated,  and  then  the  calculation  loop  of  subsection  2.8.2  is  restarted  from  the  beginning. 


l{n]  =  L[n-  1] 
L[n]  =  L 
(^Jn]  = 

t[n]  =  t[n-  1] 


(L  and  are  the  values  calculated  in  sub-subsection 


(2.85) 


2.8.3  Determination  of  specimen  failure 

See  subsection  2.5.4. 

2.9  Conclusions 


This  chapter  has  presented  procedures  for  implementing  numerical  simulations  of 
standard  experimental  tests  according  to  the  equations  of  the  CLM.  After  the  material  properties 
which  appear  in  the  model’s  equations  are  determined,  the  output  from  numerical  simulations  can 
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be  compared  to  experimental  results,  in  order  to  judge  the  model’s  value.  Simulations  can  be  used 
to  calculate  lifetimes  far  beyond  the  possible  range  of  experimental  data. 


3.  ANALYSIS  OF  COMPUTER  SIMULATIONS  FOR  SEN 
SPECIMENS  SUBJECTED  TO  CONSTANT  LOADS 

3.1  Introduction 

This  chapter  presents  the  results  and  analysis  for  numerical  simulations  of  single  edge 
notched  (SEN)  specimens  subjected  to  constant  loads.  The  driving  forces  are  broken  down  into 
sub-quantities,  and  the  behaviors  of  the  sub-quantities  are  shown.  The  simulation  input 
parameters  are  given  in  dimensionless  form,  which  reduces  the  number  of  simulation  degrees  of 
freedom,  and  allows  the  specification  of  the  criterion  for  simulation  similarity.  The  factors  which 
determine  whether  a  simulation  is  discontinuous  or  continuous  are  specified.  Simulations  using 
the  complete  Crack  Layer  Model  (CLM)  are  compared  to  simulations  using  the  simplified  CLM. 
Lifetime-stress  plots  are  shown  and  discussed  for  both  discontinuous  and  continuous  simulations. 
Lifetime-stress  plots  are  compared  for  discontinuous  simulations  of  tension,  eccentric  tension  and 
pure  bending. 

3.2  Forcing  and  resisting  forces 

From  equations  2.2(a)  and  2.3,  the  driving  force  for  crack  length  extension  (X^)  can  be 

written 

Xi=  +  (3.1) 

Equation  3.1  can  also  be  expressed 

Xi  =  D,-Ri.  (3.2) 

Di  and  R/  are  defined  as 

Rl  =  <^dry^dr+^y 
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2y  represents  the  material  function  2y(Atp^)  (see  equations  2.4,  2.49,  2.50).  Recall  that  in  this 


thesis,  the  drawing  stress  is  a  negative  quantity  (equation  2.17).  Therefore,  within  this 

thesis,  the  COD  due  to  the  drawing  stress  (6^^)  is  also  a  negative  quantity.  Dj  is  called  the  forcing 

part  of  the  driving  force  for  crack  length  extension,  and  is  always  positive.  /?;  is  called  the 

resisting  part  of  the  driving  force  for  crack  length  extension,  and  is  also  always  positive. 

Recall  equation  2.2(b), 


= 


(K^  +  K^,)(K^  +  (l  +  2y)K^^) 


(3.4) 


is  the  driving  force  for  crack  layer  extension,  and  can  be  written 


X^  =  D,-R,. 


(3.5) 


and  are  defined  as 


R,=- 


E 


(3.6) 


y  is  the  material’s  specific  energy  of  drawing,  with  domain  (0  <  y  <  «>) .  Because  is  negative, 

within  this  thesis  is  always  negative.  is  called  the  forcing  part  of  the  driving  force  for 

crack  layer  length  extension,  and  is  always  positive.  7?^  is  called  the  resisting  part  of  the  driving 

force  for  crack  layer  length  extension,  and  is  also  always  positive. 

Some  statements  will  be  written  about  equations  3. 1-3.6,  which  are  always  physically  true, 
and  can  be  shown  to  be  always  numerically  true. 
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5-2(1 +yW  (3,7) 

i:_>(l+2Y)|ii:J 

If,  for  instance,  y  is  very  large,  then  the  second  of  relations  3.7  requires  that  the  magnitude  of  5^^ 
(the  COD  at  X  =  I  due  to  the  process  zone  stress  applied  on  the  crack  faces  between  x  =  I 
and  X  =  L )  be  very  small  with  respect  to  5^  (the  COD  at  x  =  I  due  to  the  remote  load  stress 
)  applied  on  the  crack  faces  between  x  =  0  and  x  =  L).  Also  if  y  is  very  large,  then  the 
fourth  of  relations  3.7  requires  that  the  magnitude  of  (the  SIF  at  jr  =  L  due  to  the  process 
zone  stress  (o^^)  applied  on  the  crack  faces  between  x  =  I  and  =  L)  be  very  small  with 
respect  to  (the  SIF  at  j:  =  L  due  to  the  remote  load  stress  (ct^(A:) )  applied  on  the  crack  faces 

between  x  =  0  and  x  =  L).  As  Figure  6  shows,  as  y  increases,  the  equilibrium  crack  layer 
length  which  is  the  upper  bound  for  the  crack  layer  length  (L),  decreases,  and  L-l 

decreases.  As  L-l  decreases,  the  ratios  and  decrease  fast  enough  so  that  the 

second  and  fourth  of  relations  3.7  are  always  true. 

In  Figure  16,  the  crack  length  (/  )  and  the  crack  layer  length  (L  )  have  values  so  that  the 
forcing  part  of  the  driving  force  for  crack  layer  length  extension  (Dj)  equals  zero.  From  the  first 
of  equations  3.6,  this  implies  that 


KJL  ,...)  + ,L  =  0. 


(3.8) 
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In  Figure  16,  the  dimensionless  (barred)  quantities  are  normalized  by  the  values  of  the  quantities 
*  * 

using  I  and  L  ,  i.e., 


Diil,  L) 
W^ihL) 


Did,  L) 

♦  * 
Did  ,  L  ) 

DLd,  L) 

*  * 


(3.9) 


Figure  16(a)  shows  the  evolution  of  and  as  I  is  changed  from  I*  to  L*  while  L  is  held 

fixed  and  equal  to  L  .  Figure  16(b)  shows  the  evolution  of  and  as  L  is  changed  from  /* 

*  ,  ,  *  - 

to  L  while  /  is  held  fixed  and  equal  to  /  .  Figure  16(b)  also  shows  the  evolution  of  when  the 

crack  layer  length  is  close  enough  to  the  specimen  width  ( VF) ,  so  that  the  edge  influence  is 

significant.  Note  as  Figure  16  shows,  that  when  /  =  L  (in  Figure  16  when  /  =  L  or  L  =  /  ),  £)/ 

equals  zero.  When  the  crack  tip  and  the  front  of  the  crack  layer  have  the  same  location  (/  =  L) , 
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then  the  crack  tip  opening  displacements  (CTODS),  i.e.,  5^  and  6^^  each  have  value  zero.  When 
the  CTODS  are  equal  to  zero,  then  the  first  of  equations  3.3  shows  that  (or  )  equals  zero. 


Figure  16.  Forcing  forces  as  functions  of;  (a)  crack  length;  (b)  crack  layer  length 


3.3  Dimensionless  representation  of  simulation  esiput  parameters 

The  goal  of  this  section  is  to  represent  all  simulation  input  parameters  in  dimensionless 
form.  By  doing  this,  the  number  of  simulation  degrees  of  freedom  is  reduced,  and  the  criterion  for 
simulation  similarity  can  be  specified. 

The  material  constants  which  are  used  as  normalizing  factors,  are  Young’s  modulus 
(  force  ^ 

E,  dimensions  = - -  ,  the  energy  per  unit  crack  surface  area  for  drawn  but  undegraded 

^  (length)  ^ 
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material  [lyn,  dimensions  =  the  degradation  coefficient  f A:rv,  dimension  = ,  and 

V  lengthy  °  ^  timey 


the  drawing  stress  dimensions  = 


force 

(length)' 


The  system  characteristic  length  (/*)  is 


defined  as 


L  = 


IJqE 

2  ■ 
<^dr 


(3.10) 


Note  that  the  dimension  of  /*  is  length.  The  dimensionless  simulation  input  parameters  can  now 
be  specified. 

The  dimensionless  specimen  width  (W)  is 


W  = 


W 

/*’ 


(3.11) 


( W  dimension  =  length) . 


The  dimensionless  crack  length  kinetic  coefficient  (k,)  is 


k  -  —k 
1  1’ 


(3.12) 


(k^  dimensions  =  — ') 

V  ^  (force) (time) j 


The  dimensionless  crack  layer  length  kinetic  coefficient  (ik,)  is 


k  -  —k 

'^2  “  t  e’'^2’ 


kQE 


(3.13) 


dimensions  = 


(length)^  A 
(force)  (time)  j ' 
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The  dimensionless  specific  energy  of  drawing  (y)  is 

y  =  y, 

(y  dimension  =  1 ) . 

The  dimensionless  initial  crack  length  (/q)  is 


(/q  dimension  =  length) . 

The  dimensionless  initial  crack  layer  length  (Lq)  is 


(Lq  dimension  =  length) . 

The  dimensionless  initial  time  increment  (Atg)  is 

A^o  =  ^qATq, 

(AtQ  dimension  =  time) . 

The  dimensionless  remote  load  stress  (o„(jc))  is 


Ooo(^) 


<yoo(^) 


dimensions  = 


force 

(length)^ 


(3.14) 


(3.15) 


(3.16) 


(3.17) 


(3.18) 


A  particular  dimensionless  simulation  is  defined  by  specifying  the  set  of  dimensionless  input 
parameters  {W,kpk2,y,lQ,LQ,At,a^{x)}.  The  set  of  normalizing  parameters 
{E,  IJq,  kQ,  is  eliminated,  and  these  simulation  input  parameters  would  each  be  set  equal  to 
1.  For  a  dimensionless  simulation,  every  calculated  quantity  is  dimensionless,  i.e.,  stress  intensity 
factors  (SIFs),  crack  opening  displacements  (CODs),  values  of  the  energy  release  rate  (7j) ,  etc. 

Simulation  “A”  and  “B”  are  similar  if  their  sets  of  dimensionless  input  parameters  are  identical. 
For  similar  simulations,  outputs  can  be  scaled.  For  instance,  if  simulations  “A”  and  “B”  are 
similar,  then  the  total  time  for  simulation  “B”  can  be  found  from  the  total  time  for 

simulation  ‘A”  (ty(A))  >  without  having  to  aetually  perform  simulation  “B”.  If  the  dimensionless 
total  simulation  times  are  equal,  then 


V(A)  -  V(B)- 

Therefore, 


(3.19) 


V(B)  - 


^O(A) 

^O(B) 


V(A)  • 


(3.21) 


3.4  Discontinuous  and  continuous  growth  modes 

Figure  17  shows  graphs  made  from  two  computer  simulations  of  an  SEN  specimen 
subjected  to  creep  loading.  Each  graph  plots  the  time  evolutions  of  the  crack  length  (/),  crack 
layer  length  (L),  and  equilibrium  crack  layer  length  (L^^).  Figure  17(a)  shows  a  simulation  in 
which  the  growth  rates  of  the  three  plotted  variables  are  always  positive,  so  this  type  of  growth  is 
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called  continuous.  This  means  that  the  forcing  part  of  the  driving  force  for  crack  length  extension 
(D;)  is  greater  than  the  resisting  part  of  the  driving  force  for  crack  length  extension  (R[),  during 

the  entire  simulation.  In  Figure  17(b),  intervals  of  positive  growth  rates  alternate  with  intervals  of 
zero  growth  rates,  so  this  type  of  growth  is  called  discontinuous.  This  means  that  time  intervals  in 
which  D[ >  alternate  with  time  intervals  in  which  Di<Ri.  Note  that  the  time  scales  used  in 
Figures  1 1(a)  and  1 1(b)  are  not  equal. 


a)  continuous 


Figure  17.  Crack  layer  growth 


Simulation  “A”  is  identical  to  simulation  “B”,  if  for  each  simulation  input  parameter  “C”,  the 
value  of  “C”  for  simulation  “A”  equals  the  value  of  “C”  for  simulation  “B”.  Simulation  “A”  is 
comparable  to  simulation  “B”,  if  exactly  one  simulation  input  parameter  prevents  “A”  and  “B” 
from  being  identical.  If  a  particular  simulation  is  discontinuous,  then  a  continuous  comparable 
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simulation  can  be  made  by  either  increasing  the  input  parameter  corresponding  to  the  remote  load 
(a^(x)),  or  by  adjusting  the  input  parameter  associated  with  temperature  {2'^q)  to  reflect  a 

temperature  increase.  The  material  constant  2Yq  is  called  the  initial  surface  energy,  and  represents 
the  rate  of  energy  required  for  extending  the  crack  in  undegraded  drawn  material  (i.e.,  in  material 
which  has  been  drawn  for  a  time  interval  At  =  0).  is  a  decreasing  function  of  temperature. 
So,  a  comparable  discontinuous  simulation  can  be  made  from  a  continuous  simulation,  by  either 
increasing  a^(x)  or  by  decreasing  2Yo  •  If  from  the  start,  a  simulation  is  discontinuous,  it  means 


D^(t  =  0)  is  less  than  =  0)  (equations  3.2,  3.3).  If  a^(x)  is  increased,  then  £);(t  =  0)  and 


/?y(t  =  0)  will  both  increase,  but  Z)^(t  =  0)  will  increase  more.  can  be  increased  enough 

so  that,  Di{t  =  0)  is  greater  than  =  0) ,  and  the  simulation  becomes  initially  continuous.  If 
2Yo  is  decreased,  then  =  0)  will  decrease,  while  =  0)  will  not  change.  2Yo  can  be 
decreased  enough  so  that  the  simulation  would  also  become  initially  continuous. 

As  a  simulation  continues,  the  crack  length  (/)  and  the  crack  layer  length  (L)  will  grow.  As  I 
and  L  increase,  the  quantity  +  ^dr  (5^^  is  negative)  will  increase.  This  means  that  all 
simulations  (if  not  manually  stopped)  will  end  as  continuous.  Using  equations  3.3,  when 


(3.22) 


becomes  true,  the  remainder  of  the  simulation  will  be  continuous.  For  a  simulation  which  is 
entirely  continuous,  condition  3.22  will  be  true  at  time  t  =  0 . 


79 


In  Figure  17(a),  for  every  time  t,  the  crack  layer  length  (L)  is  smaller  than  the  equilibrium 
crack  layer  length  In  Figure  17(b),  for  most  times  1,  L  is  almost  equal  to  (L  can 

approach  closely  to  L  ,  but  is  always  prohibited  from  reaching  it).  The  ratio  of  the  crack  layer 
length  kinetic  coefficient  (^2)  to  the  crack  length  kinetic  coefficient  (k^)  determines  how  closely 
L  can  approach  to  L  .  The  observation  of  computer  simulations  has  shown  that  if  the  ratio 
(it2/itj )  of  these  two  material  constants  is  at  least  1/10 ,  then  it  is  possible  for  L  to  approximate 

^eq- 

For  a  continuous  simulation,  the  lifetime  (ty)  is  strongly  dependent  on  kj  and  ^2-  If 
k,  =  k2 ,  then  tjr  is  approximately  inversely  proportional  to  them.  For  a  discontinuous  simulation 
(a  simulation  in  which  condition  3.22  becomes  true  just  before  the  simulation  finishes)  fy  is  only 
weakly  dependent  on  kj  and  k2.  kj  and  k2  only  affect  the  sizes  of  the  time  intervals  when 

growth  occurs,  and  for  a  discontinuous  simulation,  these  intervals  are  usually  designed  to  be  much 
smaller  than  the  time  intervals  without  growth.  For  a  discontinuous  simulation  using  a  particular 
value  for  the  remote  load  o^{x) ,  depends  mostly  on  the  material  constants  2yQ  and  ko  (ko  is 

the  degradation  coefficient,  see  equations  2.49,  2.50).  These  constants  determine  the  sizes  of  the 
time  intervals  without  growth. 

During  the  time  periods  between  intervals  of  discontinuous  growth,  the  values  of  the  driving 
forces  for  crack  length  extension  (X^)  and  crack  layer  length  extension  (X^)  equal  zero  (see 

equation  2. 1).  During  intervals  of  discontinuous  growth,  and  are  positive.  So  for  a  cycle  of 
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no  growth  -  growth  -  no  growth ,  the  driving  force  rates  (X^  and  )  increase  from  zero,  and  then 
decrease  back  to  zero.  This  type  of  growth  is  called  stable.  When  a  discontinuous  simulation 
becomes  continuous,  and  X^  monotonically  increase.  The  growth  is  then  said  to  have  become 

unstable.  For  a  simulation  which  is  continuous  from  t  =  0  (a  continuous  simulation),  X^  and 

are  always  monotonically  increasing,  and  all  growth  is  unstable. 

Most  experimental  observations  have  been  of  discontinuous  growth,  and  have  yielded  a  linear 
relationship  between  the  logarithms  of  the  applied  stress  and  lifetime  (see  Figure  18),  i.e., 

log(r^)  =  a  +  piog(aJ.  (3.23) 

The  experimentally  observed  interval  (range)  of  coefficient  P  is  approximately  (-2.5, -5.5).  It 
will  be  shown  in  section  3.7  that  discontinuous  simulations  can  be  used  to  reproduce  this  interval, 
while  continuous  simulations  cannot. 
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Figure  18.  Experimental  stress-lifetime  relationship 


3.5  STMPT.TFIKn  MODEL  FOR  DISCONTINUOUS  GROWTH 

The  simplification  used  to  implement  the  simplified  CLM  for  discontinuous  growth,  is  that 
the  crack  layer  length  (L)  is  always  equal  to  the  equilibrium  crack  layer  length  (see  equation 
2.53).  This  simplification  results  in  growth  intervals  which  happen  in  zero  time  (see  Figure  19). 
Instantaneous  growth  intervals  could  be  instituted  within  numerical  simulations  of  the  complete 
CLM  if  it  was  possible  to  set  each  of  the  kinetic  coefficients  (/tj  and  ^2  >  see  equation  2.1)  equal 
to  infinity. 


Figure  19.  Discontinuous  growth  using  simplified  CLM 


Figure  20  is  a  plot  of  lifetime  (tj-)  versus  the  kinetic  coefficients  and  k2  (for  this  plot,  the 
degradation  coefficient  equals  1,  see  equation  2.51)  made  from  simulations  using  the  complete 
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CLM.  The  horizontal  lines  are  the  lifetimes  for  the  same  simulations  using  the  simplified  CLM. 
The  complete  and  simplified  lifetimes  become  almost  identical  when  and  k2  are 
approximately  equal  to  10.  Experimental  observations  have  shown  growth  corresponding  to  large 
values  of  and  k2.  So  the  simplified  model  can  be  used  to  accurately  predict  the  lifetimes  of 

real  specimens.  All  calculations  for  plots  which  appear  in  the  rest  of  this  chapter  are  made  from 
simulations  using  the  simplified  CLM. 


0  2  4  6  8  10 
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Figure  20.  Lifetime  versus  kinetic  coefficients 


Figure  21  is  a  plot  made  from  computer  simulations  of  the  dimensionless  critical  (failure) 
crack  length  (/^)  versus  the  logarithm  of  the  dimensionless  specimen  width  (IT) .  (T^)  and  (W) 
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are  normalized  by  /*  (see  equation  3.10).  The  values  of  (/^)  at  the  right  of  the  plot  are  almost 
identical  to  those  for  the  analytical  solutions  for  a  semi-infinite  plate. 


log.o  (W) 


Figure  21.  Failure  crack  length  versus  specimen  width 


3.6  Lifetimes  for  tension  simulations 

Figures  22(a)  and  (b)  are  plots  of  the  logarithm  of  dimensionless  lifetime  {t^)  versus  the 

logarithm  of  applied  remote  stress,  made  from  simulations  of  specimens  subjected  to  constant 
tension  loads.  Figure  22(a)  is  made  from  simulations  of  discontinuous  growth,  and  Figure  22(b)  is 
made  from  simulations  of  continuous  growth.  Best  fit  lines  are  drawn  thru  the  points  for  each  of 

the  values  of  the  dimensionless  specimen  width  ( W) .  The  only  simulation  input  parameter  which 
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Figure  22.  Tension  simulations  lifetime-stress  relations 


varies  in  Figures  22  is  W .  Since  Iq/W  is  constant,  as  W  changes  the  specimen  is  geometrically 

scaled.  The  formal  explanation  for  why  the  dimensionless  lifetime  (tj.)  changes  as  the  specimen 
changes  size  is  that  the  similarity  referred  to  in  section  3.3  is  violated.  The  physical  explanation  is 
that  geometrically  scaling  a  specimen  also  scales  the  specimen’s  energy  release  rate  (J\).  For 

instance,  if  the  size  of  a  specimen  is  doubled,  then  so  is  .  But  changing  the  size  of  a  specimen 
does  not  change  the  properties  of  the  material  from  which  the  specimen  is  made,  i.e.,  the 
specimen’s  resistance  to  crack  extension  (the  material  function  is  constant  for 

geometric  scaling.  A  crack  will  grow  when  the  driving  force  for  crack  length  extension 
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{Xi  =  7i  -2j(Atp^))  becomes  positive.  As  W  increases,  increases,  the  waiting  time  until 

each  interval  of  crack  growth  decreases,  and  the  lifetime  (r^)  decreases.  W  =  32  is  used  in 

Figure  22(b),  because  32  is  very  close  to  being  the  smallest  value  of  W  for  which  the  simulation 
is  initially  continuous  (remember  that  within  this  thesis,  a  continuous  simulation  is  one  which  is 
continuous  at  time  t  =  0,  and  a  discontinuous  simulation  is  one  which  is  discontinuous  at  time 

t  =  0).  So  W  =  32  is  approximately  the  lower  limit  of  W  for  a  continuous  simulation.  Note  in 
Figure  22(b)  that  as  W  is  increases,  the  best  fit  lines  appear  to  converge  to  a  limit. 
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Figures  23  are  plots  of  the  coefficients  a  and  |3  (see  equation  3.23)  of  the  best  fit  lines  of 
Figures  22.  Note  that  the  discontinuous  simulations  can  reproduce  the  experimentally  observed 
interval  for  P ,  (-2.5,  -5.5) ,  while  the  continuous  simulations  cannot. 


>og,o(M^) 

a)  Discontinuous  simulations 


Figure  23.  Log(lifetime)  -  Log(applied  stress)  best  fit  line  coefficients 
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3.7  I.IFETIMES  FOR  ECCENTRIC  TENSION  AND  PURE  BENDING 

This  section  compares  lifetimes  for  discontinuous  simulations  of  tension,  eccentric  tension 
and  pure  bending.  The  input  parameters  are  adjusted  so  that  each  type  of  loading  has  the  same 
maximum  applied  stress  (see  Figure  24). 


e~— 


Figure  24.  Applied  stress  distributions 


Figures  25  show  the  calculated  lifetime-applied  stress  relationships.  Note  that  for  the  small 
value  of  W  the  coefficients  P  of  equation  3.23  are  approximately  equal,  while  for  the  big  value  of 
W  each  value  of  P  is  distinct.  (Note  that  the  “tension”  best  fit  line  of  Figure  25(a)  and  the  top  best 
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fit  line  of  Figure  22(a)  are  the  same.  Also  note  that  the  “tension”  best  fit  line  of  Figure  25(b)  and 
the  bottom  best  fit  line  of  Figure  22(a)  are  the  same.) 


log.o('/)  10g,o(f/) 


Figure  25.  Lifetime-stress  relations  for  tension,  eccentric  tension  and  pure  bending 


3.8  Conclusions 

This  chapter  has  analyzed  the  output  from  numerical  simulations  of  SEN  specimens 
subjected  to  constant  loads.  The  behaviors  of  the  driving  forces  were  described.  The  simulation 
input  parameters  were  presented  in  dimensionless  form.  Discontinuous  and  continuous 
simulations  were  examined  and  compared.  Simulations  using  the  complete  CLM  were  compared 
to  simulations  using  the  simplified  CLM.  Simulation  lifetime-stress  relations  were  shown  to  have 
the  same  form  as  those  obtained  experimentally.  It  was  shown  that  discontinuous  simulations  can 
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reproduce  the  same  range  of  coefficient  values  in  the  lifetime-stress  relations  as  have  been 
experimentally  determined.  It  was  shown  that  continuous  simulations  cannot  reproduce  this 
range.  Lifetime-stress  relations  were  compared  for  discontinuous  simulations  of  tension,  eccentric 


tension  and  pure  bending. 


4.  RELIABILITY  CALCULATIONS  FOR  PLASTIC 
STRUCTURAL  COMPONENTS 

4.1  Introduction 

This  chapter  will  show  a  method  for  calculating  the  reliability  of  statically  loaded  plastic 
structural  components.  For  a  particular  stress  (ct)  and  time  (t),  the  reliability  is  the 

probabihty  that  the  component  is  still  functioning,  i.e.,  that  it  has  not  failed.  The  use  of  plastic  in 
engineering  structures  is  always  increasing,  so  with  respect  to  safety  and  liability  it  is  very 
important  to  be  able  to  know  the  expected  lifetime  for  a  plastic  part.  Classical  theories  have  been 
developed  for  elastic,  brittle  materials  which  can  be  used  to  calculate  reductions  in  ultimate 
strength  when  flaws  are  present.  But  recent  studies  of  plastic  field  failures  show  that  a  plastic  part 
will  develop  a  crack  after  being  loaded  for  a  long  time  by  a  relatively  small  stress  [46,47].  The 
crack  originates  from  the  part’s  largest  defect  (assuming  homogeneous  stress).  The  defect  locally 
magnifies  the  stress,  a  process  zone  develops,  and  slow  crack  growth  begins.  The  defects  are 
introduced  into  the  plastic  material  during  the  manufacturing  process,  and  are  in  the  form  of  voids 
or  dust  particles.  Tests  are  performed  on  sample  specimens  to  determine  the  probability  density 
function  (PDF)  of  the  critical  (largest)  defect  size.  The  defect  is  considered  to  be  a  crack  with 
length  equal  to  the  defect  size.  For  a  particular  applied  stress  and  structural  component  geometry, 
the  equations  of  the  crack  layer  model  provide  a  relation  between  the  initial  crack  length  (defect 
size)  and  lifetime.  Having  both  the  PDF  for  initial  crack  length  (defect  size)  and  the  equation 
linking  initial  crack  length  and  lifetime,  allows  the  calculation  of  the  PDF  for  lifetime.  Once  the 
PDF  for  lifetime  is  known,  the  reliability  function  for  the  structural  part  can  be  determined. 

4.2  Method  assumptions 

1)  The  sizes  and  locations  of  the  defects  are  random  and  cause  local  stress  concentrations. 
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2)  The  stress  concentration  induced  by  a  defect  results  in  process  zone  formation. 

3)  A  crack  in  a  structural  part  originates  from  the  largest  defect.  Therefore  one  of  the  PDFs 
from  the  distribution  of  extremes  should  be  used  for  the  PDF  of  defect  size. 

4)  The  crack  layer  model  can  adequately  predict  lifetimes  for  plastic  specimens. 

4.3  Analysis  of  defects 

Destructive  tension  tests  in  which  the  load  is  increased  until  failure  occurs  are  performed 
on  many  smooth  (unnotched)  specimens  of  the  given  material.  The  fracture  surface  of  each 
specimen  is  examined,  and  the  size  of  the  defect  from  which  the  crack  originated  (which  is 
assumed  to  be  the  largest  defect  in  the  specimen)  is  measured.  A  histogram  is  constructed  from 
the  failure  defect  sizes,  and  the  parameters  of  the  PDF  of  the  statistics  of  extremes  for  a  finite 
interval  are  fitted  to  the  histogram.  The  cumulative  distribution  function  (CDF)  of  the  statistics  of 
extremes  for  a  finite  interval  is  [48] 


FHo)  = 


^  *  ^0  ^  ^min 


1  :  L..  <  I 


-max 


‘0 


(4.1) 


The  PDF  (/(/q))  is  obtained  by  differentiating  equation  4.1  with  respect  to  Iq  (the  initial  crack 

length  or  defect  size).  Figure  26  shows  /(/q)  for  an  observed  set  of  critical  defect  sizes,  Iq  ,  made 

from  the  histogram  of  destructive  tension  tests  performed  on  polycarbonate  bars,  and 

(equation  4.1)  are  the  minimum  and  maximum  values  from  the  histogram,  a  and  P  are  chosen  so 
that  the  PDF  most  closely  approximates  the  shape  of  the  histogram.  Their  are  various  statistical 
criteria  which  can  be  used  to  judge  the  best  fit. 
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Figure  26.  Modeling  probability  density  function  for  maximum  defect  size 


4.4  KP4ETICS  OF  CRACK  LAYER  EVOT.UTTON 

The  driving  force  for  crack  extension,  ,  is  the  derivative  of  Gibb’s  potential,  G ,  with 

respect  to  the  crack  length,  / .  The  driving  force  for  crack  layer  length  extension,  ,  is  the 

derivative  of  G  with  respect  to  the  crack  layer  length,  L  [40].  For  slow  (quasi-static)  crack 
growth,  the  rates  of  growth  of  the  erack  and  crack  layer  are  assumed  to  be  linearly  proportional  to 
the  driving  forces. 
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L  =  k.X^ 

Evaluation  of  the  driving  forces  requires  knowledge  of  material  properties,  such  as  the  drawing 
stress,  draw  ratio,  specific  energy  of  drawing.  Young’s  modulus  and  the  characteristic  time  of 
material  degradation  within  the  process  zone.  Also  required  are  conventional  fracture  mechanics 
computations  of  the  stress  intensity  factors  and  crack  opening  displacements. 

The  system  of  equation  4.2  consists  of  two  strongly  bound  non-linear  ordinary  differential 
equations,  and  can  only  be  solved  numerically.  Numerical  simulation  of  the  model  at  moderate 
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stress  levels  results  in  discontinuous  crack  growth,  similar  to  experimentally  observed  behavior 
(Figure  27). 
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Figure  27.  Discontinuous  growth  of  crack  and  process  zone 


4.5  Lifetime  -  initial  defect  .sizk  rft.atiqnship 

Computations  of  the  total  time  until  part  failure,  (lifetime),  produce  a  linear 
relationship  with  initial  defect  size,  Iq 

V  =  =  A  +  51n(/o).  (4.3) 

The  coefficients  (A,  5)  of  equation  4.3  depend  on  the  specimen  geometry  and  type  of  loading, 
i.e.,  uniform  tension,  three  point  bending,  etc.  Figure  28  is  a  calculated  plot  at  various  stress  levels 
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of  lifetime  as  a  function  of  initial  defect  size  for  single-edge  notched  specimens  of  polyethylene 


piping  material.  (Note  that  a  =  — .) 

^dr 


-1  -0.5  0  0.5 

logio(^o)  (log, /mm.)) 

Figure  28.  Lifetime  versus  logarithm  of  initial  crack  length  for  five  stress  levels 


4.6  Lifetime  probability  density  function 

Since  the  lifetime  is  a  function  of  the  random  initial  defect  size  (/q)  (equation  4.3), 

and  since  the  PDF  for  initial  defect  size  (/(/q)  )  (derivative  of  equation  4.1)  is  known,  a  standard 
transformation  from  probability  theory  [49]  allows  the  calculation  of  the  PDF  for  lifetime 
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gitf)  = 


d(tf) 


(4.4) 


Figure  29  shows  the  PDF  for  lifetime  for  various  normalized  stress  levels.  (Figure  29  uses  the  data 
from  Figures  26  and  28.  This  means  that  Figure  29  shows  PDFs  for  a  hypothetical  material.) 


0.1  1  10 
TIME  (years) 

Figure  29.  Lifetime  probability  density  functions  for  five  stress  levels 


4.7  Reliability  of  structural  component 

The  reliability  function,  which  gives  the  probability  that  a  part  is  still  functioning  at  a 
given  time  t,  depends  on  the  stress  level  and  is  defined  as 
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R{t,G)  =  1-  J  8a(tf)ditf).  ■  (4.5) 

V  =  o 

Figure  30  is  a  plot  of  lifetime  versus  remote  stress  for  various  values  of  reliability.  (Figure  30  uses 
the  data  of  Figure  29.) 


Figure  30.  Lifetime  versus  stress  level  for  four  values  of  reliability 
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A  simpler  measure,  the  mean  time  to  failure,  (tjr) ,  is  shown  in  Figure  31  versus  stress  level. 
(Figure  31  uses  the  data  of  Figure  29.) 


Figure  3 1 .  Mean  lifetimes  for  five  stress  levels 


4.8  Conclusions 

A  method  of  reliability  detemaination  for  plastic  structural  components  has  been 
presented.  The  method  is  based  on  the  statistical  determination  of  the  probability  density  function 
for  initial  defect  size,  and  on  the  connection  provided  by  the  crack  layer  model  between  initial 
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defect  size  and  lifetime.  The  method  can  be  used  to  recommend  replacement  intervals  for 
structural  parts,  assist  in  material  selection,  and  decide  on  design  strategies. 


5.  COMPUTER  SIMULATION  OF  PLASTIC 
TOUGHNESS  TEST 

5.1  Introduction 

Since  plastics  are  increasely  being  used  instead  of  metals  for  industrial  applications,  there 
is  an  interest  for  obtaining  a  standard  method  of  quantifying  the  fracture  toughness  of  a  plastic. 
The  first  criterions  used  were  simple  empirical  tests  which  gave  little  useful  information.  More 
recently,  efforts  have  been  made  to  apply  metal  testing  procedures  to  plastics.  An  ASTM 
subcommittee  (D-20)  was  created  to  develop  methods  to  rate  plastics  based  on  resistance  to  crack 
initiation  (or  fracture  toughness)  [50],  (ASTM  Standards  E8 13-89,  El  152-89).  First,  there  were 
attempts  to  find  the  linear  elastic  fracture  parameters  K^c  ^tid  (which  are  used  as  measures 

of  metal  fracture  toughness).  Round-robin  tests  were  performed  at  several  laboratories.  Because 
plastic  response  is  usually  nonlinear  and  inelastic,  these  parameters  were  not  determinable. 
Attempts  were  then  made  to  find  the  elastic-plastic  fracture  entities  and  the  7j  curve.  Round- 

robin  tests  were  again  conducted,  and  it  was  found  that  these  measures  were  not  able  to  be 
uniquely  determined,  and  so  cannot  be  material  properties. 

In  this  chapter,  the  crack  layer  model  will  be  used  to  find  and  the  7j  curve  as  specified  in 
ASTM  E8 13-81.  It  will  be  shown  that  these  measures  are  not  unique.  Then,  the  (material)  time 
until  failure  (or  lifetime)  ty ,  as  found  from  the  solution  of  the  equations  of  the  crack  layer  model 

will  be  proposed  as  a  measure  of  fracture  toughness. 

5.2  The  crack  layer  model 

The  crack  layer  model  (CLM)  was  derived  from  the  observation  that  for  many  engineering 
plastics,  the  material  near  the  path  of  a  slowly  growing  crack  is  permanently  altered  when  the 
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crack  tip  passes  by  [43].  The  transformed  material  is  called  the  process  zone  (PZ).  The  change  in 
the  PZ  material  is  that  it  becomes  drawn  (permanently  stretched)  in  the  direction  of  the  applied 
loading  (perpendicular  to  the  crack).  The  drawing  happens  because  the  crack  tip  magnifies  the 
stress  in  the  material  nearby,  and  the  (limiting)  drawing  stress  (cr^^)  is  reached  within  this 
material.  Figure  32  shows  using  a  single  edge  notched  (SEN)  specimen,  that  the  crack  length  is 
called  I ,  and  the  process  zone  length  is  called  .  The  crack  and  process  zone  are  together  called 
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the  crack  layer.  The  length  of  the  crack  layer  is  named  L .  For  a  material  capable  of  exhibiting  a 
process  zone,  the  process  zone  will  always  exist  ahead  of  the  crack  tip  (assuming  a  nonzero  stress 


102 


field).  As  the  crack  grows,  is  not  constant.  Therefore,  I  and  L  are  not  always  equal.  Because  of 

this  it  seems  natural  to  treat  I  and  L  as  independent  variables.  This  is  what  is  done  in  the  CLM.  In 
this  model,  the  two  independent  variables  are  /(r,  L,  . . . )  and  ...),  and  each  is  a  function  of 
the  other.  The  driving  forces  for  crack  length  extension  (Z^)  and  crack  layer  length  extension 
(Z^)are 

X,  =  7,-2Y(A/p,)  (5.1) 

(if  <  0,  then  Z;  is  set  to  0), 

+  ( 1  +  2y)X*)  (5.2) 

(if  <  0,  then  Z^  is  set  to  0). 

The  instantaneous  growth  rates  for  the  crack  and  crack  layer  are 
i=k,X, 

.  •  (5.3) 

L  =  k2Z^ 

Equations  5.3  are  a  system  of  two  nonlinear  ordinary  differential  equations,  which  can  be 
solved  simultaneously  thru  time,  to  provide  the  chronological  evolution  of  I  and  L.  For  a 
computer  simulation,  the  material  properties,  specimen  geometry,  initial  crack  length  (/g),  test 

type  and  loading  conditions  are  specified.  The  developments  of  /  and  L  are  recorded,  until  L 
reaches  the  specimen  width  ( W)  indicating  specimen  failure.  It  has  been  shown  that  the  numerical 
implementation  of  the  CLM  can  imitate  the  discontinuous  growth  which  has  been  observed 
during  experiments  on  many  plastics  [45]. 
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5.3  Simulation  of  displacement  control  test 

Figure  33  was  generated  from  the  numerical  simulation  of  a  displacement  control  test 

using  the  equations  of  the  CLM.  It  shows  a  typical  plot  of  remote  stress  (o^)  (and  the  associated 
plots  of  /  and  L)  versus  displacement  (A).  The  specimen  stiffness  (the  slope  of  the  vs.  A 

curve)  is  a  decreasing  function  of  I ,  and  so  has  its  maximum  value  at  the  start  of  the  simulation, 
when  /  =  /q  and  t  =  0.  The  stiffness  decreases  until  it  is  almost  zero  when  L  reaches  the 

specimen  width  W,  t  =  tj,  and  the  simulation  ends.  The  growth  of  /  is  intermittent  and 

discontinuous.  The  sizes  and  frequencies  of  growth  jumps  in  I  are  proportional  to  the  current 
value  of  the  energy  release  rate  (7j ).  For  a  displacement  control  test,  7j  has  initial  value  zero, 

reaches  its  maximum  value  at  some  intermediate  value  of  A ,  and  then  decreases.  At  each  instant 
when  /  increases,  the  equations  of  the  CLM  require  that  decreases,  resulting  in  the 


discontinuities  of  .  The  ratio  of  the  decrease  of  at  a  crack  jump,  to  the  increase  of 
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between  crack  jumps,  is  an  increasing  function  of  t  (or  equivalently  of  A).  attains  its 
maximum  value  when  this  ratio  equals  one. 


Figure  33.  Remote  stress,  crack  length  and  crack  layer  length  vs.  displacement 


5.4  ASTM  PROCEDURE 


According  to  ASTM  E8 13-81,  J  j  can  be  calculated  for  any  time  t  during  an  experiment 


from  expression 
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Jiit) 


2U(t) 


(5.4) 


In  equation  5.4,  U{t)  is  the  potential  energy  of  the  specimen,  and  B  is  the  specimen  thickness. 
U(t)  is  equal  to  the  area  under  the  vs.  A  curve  at  time  t,  and  l{t)  is  the  crack  length  at  time 
t.  The  7j  curve  for  an  experiment  is  plotted  versus  crack  extension  (5/).  is  found  as  the 
value  of  corresponding  to  the  intersection  of  the  curve  with  the  blunting  line.  The  blunting  line 
is  given  by  7 j  =  2a^5/ .  Since  for  plastics  the  yield  stress  )  cannot  be  measured  precisely,  the 

drawing  stress  can  be  substituted  for  it  in  the  blunting  line  expression. 

Figure  34  was  generated  from  numerical  simulations  of  displacement  control  tests  using  the 
equations  of  the  CLM.  It  shows  a  partial  plot  of  vs.  A  for  each  of  three  values  of  A .  For  a 

particular  value  of  A ,  increasing  A ,  decreases  l .  Decreasing  t ,  increases  the  average  value  of  2y , 
the  resistance  to  crack  extension.  Increasing  the  average  value  of  2y ,  decreases  I .  Decreasing  / , 
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increases  the  specimen  stiffness,  ct^/A.  Increasing  o^/A,  increases  a„.  Therefore,  for  a 


particular  value  of  A ,  increasing  A ,  increases  . 


Figure  34.  Remote  stress  vs.  displacement  for  three  values  of  displacement  rate 


107 


Figure  35  shows  the  7j  vs.  6/  (or  J-R)  curve,  calculated  using  ASTM  E813-81,  for  each  of 
the  vs.  A  curves  of  Figure  33.  The  blunting  line  is  shown,  and  each  of  the  three  values  of 
is  indicated. 


Figure  35.  Energy  release  rates  for  three  values  of  displacement  rate 


5.5  Conclusions 

The  J-R  curve  and  J  j  ^  are  rate  dependent,  and  so  cannot  be  material  properties.  They  do 

not  provide  a  measure  of  material  fracture  toughness.  No  universal  measure  of  material  fracture 
toughness  exists.  Toughness  depends  on  the  material  properties,  along  with  the  specific  conditions 
of  use,  i.e.,  geometry  and,  type  and  magnitude  of  loading.  One  material  can  be  said  to  be  better 
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than  another  only  with  respect  to  a  particular  set  of  usage  conditions.  The  CLM  can  be  used  to 
simulate  the  behavior  of  a  plastic  specimen  for  a  particular  set  of  usage  parameters.  Using  the 
CLM,  the  lifetimes  (tj)  for  two  materials  subject  to  identical  conditions  can  be  compared,  and 
one  material  judged  superior. 


CONCLUSIONS  AND  FUTURE  RESEARCH 


The  introduction  described  the  phenomenon  of  slow  crack  growth  in  plastics,  and  the  methods 
which  have  been  used  to  model  it. 

Chapter  1  gave  a  basic  description  of  the  theory  of  the  CLM. 

Chapter  2  presented  numerical  algorithms  for  modeling  slow  crack  growth  in  plastics  for 
various  loading  conditions,  using  the  equations  of  the  CLM.  Simplified  algorithms  were  also 
presented,  which  give  accurate  results  for  discontinuous  crack  growth  processes. 

Chapter  3  analyzed  the  results  of  computer  simulations  made  using  the  algorithms  of  Chapter 
2.  The  driving  forces  were  decomposed  and  their  behaviors  were  shown.  The  method  for  making 
a  dimensionless  simulation  was  given.  The  characteristics  of  the  two  types  of  crack  growth, 
discontinuous  and  continuous,  were  described  and  compared.  The  accuracy  of  simulations  using 
the  simplified  algorithms  was  examined.  Lifetime-applied  stress  relationships  were  presented  for 
both  discontinuous  and  continuous  simulations.  Lifetime-applied  stress  relationships  were 
compared  for  three  applied  stress  functions  (three  types  of  loading). 

Chapters  4  and  5  presented  applications  made  from  computer  simulations  using  the 
algorithms  of  Chapter  2.  Chapter  4  showed  a  method  for  calculating  the  reliability  function  for  a 
plastic  structural  component.  Chapter  5  gave  a  computer  imitation  of  the  ASTM  experimental 
method  for  measuring  plastic  fracture  toughness.  It  was  shown  that  the  computer  results  were 
displacement  rate  dependent.  It  was  suggested  that  the  experimental  results  would  also  be  rate 
dependent,  and  that  fracture  toughness  depends  on  many  usage  parameters,  and  so  is  not  a 
material  property. 

The  thesis  numerically  implemented  the  CLM  and  presented  the  results  of  computer 
simulations  of  the  implementation.  The  simulations  displayed  both  kinds  of  experimentally 
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observed  slow  crack  growth,  discontinuous  and  continuous.  The  simulations  gave  a  relationship 
between  lifetime  and  stress,  i.e.,  log(t^)  =  a  +  p  log(a^),  which  has  been  experimentally 
observed.  The  simulations  were  also  able  to  duplicate  the  experimentally  determined  range  of 
coefficient  p .  It  is  proposed  that  the  thesis  has  demonstrated  that  a  computer  program  which 
models  the  equations  of  the  CLM,  is  able  to  provide  realistic  predictions  of  the  evolutions  of  the 
crack  and  crack  layer,  and  of  long-term  lifetime  for  plastics  functioning  under  small  stresses.  It  is 
also  proposed  that  the  thesis  has  shown  that  various  load  types  and  histories  can  be  easily 
accommodated.  It  is  further  proposed  that  the  ability  to  numerically  imitate  the  slow  crack  growth 
in  plastics  has  many  practical  uses,  two  of  which  were  presented  in  Chapters  4  and  5. 

It  was  stated  in  the  introduction  that  viscoelastic  deformation  does  influence  lifetime  when 
displacement  is  prescribed.  Future  research  includes  accounting  for  viscoelastic  effects. 
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APPENDIX 

(algorithms  for  SIFs  and  CODs) 


Introduction 


This  appendix  gives  algorithms  for  calculating  the  stress  intensity  factors  and  crack 
opening  displacements  for  a  single  edge  notched  specimen.  (The  unit  dipole  SIF  green’s  function 
(G)  is  for  an  SEN  specimen  with  a  unit  dipole  force  applied  to  each  crack  face  and  normal  to  the 
CLC.  (The  CODs  used  in  this  thesis  (5„  and  6^^)  are  calculated  by  using  =  /  in 

equations  A.9  and  A.  10.) 

Abbreviations 


CLC  = 
COD  = 
SEN  = 
SIF  = 

Symbols 

/ 

L 

W 

Oix)  = 
= 


^dr  = 


crack  layer  centerline 
crack  opening  displacement 
single  edge  notched 
stress  intensity  factor 


crack  length,  also  x-coordinate  of  crack  tip 

crack  layer  length,  also  maximum  x-coordinate  of  process  zone 

length  of  SEN  specimen  in  direction  parallel  to  CLC 

the  stress  distribution  which  is  applied  normal  to  the  crack  faces 

Stress  distribution  which  would  be  present  along  location  of  CLC,  if  crack  layer 
did  not  exist.  Is  applied  normal  to  CLC  on  crack  faces  between  x-coordinates  0 
and  / ;  and  on  crack  imaginary  faces  between  x-coordinates  I  and  L.  (a „  could 
possibly  be  a  function  of  position  along  the  CLC,  i.e.,  if  for  instance,  the  specimen 
was  loaded  in  bending)  can  be  tensile  and/or  zero  and/or  compressive) 

material  drawing  stress  (applied  normal  to  CLC  on  crack  imaginary  faces  between 
x-coordinates  I  and  L)  (CJ^^  is  compressive  on  crack) 

x  axis  located  along  CLC,  or  x-coordinate  measured  along  this  x  axis  (notched 
edge  of  SEN  specimen  has  x-coordinate  zero) 


1.  Tada,  Hiroshi;  The  Stress  Analysis  of  Cracks  Handbook:  p.  2-27;  St.  Louis,  Missouri;  Del  Research  Cor¬ 
poration;  1974. 
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APPENDIX  (continued) 

Xj  =  x-coordinate  of  dipole  forces 

XcoD  -  x-coordinate  at  which  COD  is  desired  (0  <  x^qq  <  L) 

=  SIF  at  L  due  to 

=  SIF  at  L  due  to 

=  COD  at  X  due  to  .  Is  measured  normal  to  CLC.  For  x  between  0  and  I  ,  is 
equal  to  the  distance  between  crack  faces  at  location  x .  For  x  between  /  and  L ,  is 
equal  to  the  distance  between  crack  imaginary  faces  at  location  x . 

=  COD  at  X  due  to  Is  measured  normal  to  CLC.  For  x  between  0  and  I  ,  is 
equal  to  the  distance  between  crack  faces  at  location  x .  For  x  between  I  and  L ,  is 
equal  to  the  distance  between  crack  imaginary  faces  at  location  x . 

Stress  intensity  factor  equations 

jc  =  L 

j  Fix,L,W,cJx))dx  (A.l) 

x  =  0 

X=:L 

=  j  Fix,L,W,Gj^)dx  (A.l) 

x  =  l 

F(xj,  L,  W,  G)  =  gG(x^,  L,  W)  (A.3) 


G(xj,L,W)  =  -^(T^-T^  +  T^) 
JkL 


^  3.52(1 


4.35-5.28x^^^ 


(A.4) 


(A.5) 


(A.6) 
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APPENDIX  (continued) 


( 


(M)  3/2  ^ 

1.3-0.3(4"0  (;v) 

-  ^  +  0.83  -  1.764  ^ 


I  M\  2 


(N)  _  Xj  AN)  _  L 

■  L  ’  ^  -W 


Crack  opening  displacement  equations 


X  =  L 

^ooixcoD’  C7oo(^))  =  I  X,  E,  L,  W)dx 

x  =  0 


X  =  L 

^dMcoD^  E,  I,  L,  w,  J  H(xcod^  x,  E,  L,  W)dx 

X^l 


x  =  L 

H{xcod’  ^d’  E,L,W)  =  ^  J  G(xcoz),  W)G(xj,  x,  W)dx 


(A.7) 


(A.8) 


(A.9) 


(A.10) 


(A.11) 


a  -  vnz.\{xQQQ,  x^) 


(A.12) 


APPENDIX  (continued) 


Figure  36.  Crack  layer  diagram 


